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Chapter  I 
INTRODUCTION 

The  calculation  of  natural  frequencies  (poles)  and  natur¬ 
al  modes  (free  oscillations)  of  structures  is  a  fundamental 
problem  of  many  disciplines.  Until  recently,  the  mathemati¬ 
cal  study  of  these  parameters  has  unfortunately  been  limited 
to  canonic  geometries  which  lend  themselves  to  eigensoluticn 
by  separation  of  variable  techniques.  The  singularity  ex¬ 
pansion  method  (SEM)  removes  this  geometrical  restriction  by 
enabling  one  to  obtain  the  natural  frequencies  and  natural 
modes  of  an  arbitrary  object.  The  SEM  also  enables  one  to 
determine  the  response  of  the  object  to  an  arbitrarv  forcing 
function  directly  from  an  appropriate  expansion  of  the  modal 
and  pole  structure. 

The  basic  theoretical  foundations  of  the  SEM  were  ini¬ 
tially  presented  using  frequency-domain  techniques  applied 
to  electromagnetic  equations  by  3aum  [ 1 J .  Baum's  develop¬ 
ment  was  subsequently  extended  by  Marian  and  Latham  [2];  and 
rigorous  mathematical  justification  of  some  of  the  basic 
foundations  has  recently  been  presented  by  Earr.m  [31.  Ana¬ 
lytic  f requer.oy-domain  SEM  results  for  the  perfect  conduct¬ 
ing  spherical  scatterer  were  originally  obtained  by  Baum 
[l1,  and  numerical  frequency-domain  results  for  thin,  per- 


feet  conducting  cylindrical  surfaces  were  initially  present¬ 
ed  by  Tesche  [ 4 ] . 

Interest  in  time-domain  techniques  in  the  5SK  has  not 
been  as  widespread  as  f requency-donain  methods;  however, 
several  varied  contributions  have  recently  been  made  toward 
establishing  the  versatility  of  time-domain  methods.  A  time- 
domain  method  analogous  to  the  original  frequency- domain 

method  may  be  found  in  Baum  [5].  The  applicability  of 
this  method,  however,  has  been  somewhat  limited  due  to  the 
level  of  difficulty  of  the  describing  equations.  Van  Blari- 
cum  and  Mittra  [5!  developed  a  rather  unique  method  whereby 
the  natural  responses  may  be  obtained  using  Prcny’s  method 
’"l  once  the  transient  resocr.se  of  the  object  is  known.  An 
obvious  complication  with  this  method  is  that  the  determina¬ 
tion  of  the  transient  response  can  be  a  nen-triviai  problem. 
An  alternate  time-domain  method  which  sidesteps  the  compli¬ 
cations  of  the  above  methods  has  been  introduced  by  Cordaro 
and  Davis  [3].  This  method  ,  known  as  time-domain  SEP! 
(TD-SEM),  enables  one  to  find  the  natural  responses  directly 
from  the  finite ■ difference  representation  of  the  governing 
integral  equations  cast  in  a  matrix  eigenvalue  form.  Unfor- 
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results  which  have  been  obtained,  however,  indicate  that  the 
Cordaro-Davi s  method  is  capable  of  producing  a  great  deal  of 
information  quite  efficiently.  The  intent  of  this  study  is 
to  extend  the  applicability  of  TD-SEM,  and  extend  numerical 
time-domain  techniques  in  general. 

Toward  establishing  this  intent,  the  following  (prin'-'- 
pal)  set  of  tasks  are  defined:  (i)  determine  stability 
criteria  for  various  finite  difference  representations  . 
electromagnetic  equations,  (2)  develop  simple  time-dcma 
expressions  for  determining  the  SEM  coupling  coefficients 
(these  are  parameters  which  couple  the  natural  frequencies 
and  modes  to  the  incident  forcing  function),  (3)  develop  an 
eigensoluticn  algorithm  which  will  solve  the  large  scale  ma¬ 
trices  generated  by  TD-SEM,  (4)  obtain  a  pole  distribution 
for  the  linear  scatterer  discretized  with  a  large  number  of 
unknowns,  (5)  apply  TD-SEM  to  the  two-dimensional  rectangu¬ 
lar  plate  problem.  The  outline  for  establishing  these  tasks 
is  as  follows. 

The  fundamental  governing  equations  of  electromagnetics 
are  developed  in  integral  form  in  Chapter  2  using  dyadic 
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Chapter  II 


FUNDAMENTAL  INTEGRAL  EQUATIONS  OF 
ELECTROMAGNET ICS 


2 . 1  INTRODUCTION 

Singular  integral  equations  (or  singular  integro-dif f er- 
ential  equations)  represent  a  powerful  and  widely  used  ap¬ 
proach  to  the  solution  of  both  antenna  and  electromagnetic 
scattering  problems.  A  variety  of  methods  may  be  used  to 
obtain  these  equations.  Poggio  and  Miller  [9]  rigorously 
develop  the  necessary  results  using  the  vector  Green's  theo¬ 
rem  [10].  In  this  formalism,  the  concept  of  incident  and 
scattered  fields  ir.  conjunction  with  equivalent  sources  de¬ 
velops  in  a  natural  way.  In  this  chapter,  the  frequency-do¬ 
main  equations  are  developed  from  linear  system  foundations. 
Although,  perhaps,  this  approach  is  less  rigorous  than  the 
method  of  Poggio  and  Miller,  it  yields  fundamental  results 
readily,  without  extensive  vector  manipulations.  The  time- 
domain  representations  of  these  equations  are  then  obtained 
by  inverse  Fourier  transform  techniques.  These  general  fre- 
cruencv-  and  time-  domain  results  are  final Iv  specialised  to 
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2.2  MATHEMATICAL  FORMALISM 

The  mathematical  formulation  of  electromagnetic  phenomena 
is  fundamentally  dependent  on  a  concise  set  of  equations 
known  as  Maxwell's  equations.  The  complexity  of  these  equa¬ 
tions  is  highly  dependent  on  the  host  medium.  We  will  res¬ 
trict  our  discussion  throughout  to  a  homogeneous,  linear, 
and  isotropic  medium.  For  such  a  medium.  Maxwell's  equa¬ 
tions  may  be  written  in  differential  form  in  the  frequency- 
domain  as  (a  vector  will  be  denoted  by  a  single  bar;  a  fre- 
quency-dcmain  quantity  will  be  denoted  by  a  elide) 

7  ;c  ( r ;  jj )  =  -juu  Hu  ( r ;  u>) 

7  x  liC  (r;a)  =  juj£Q  EC  (r;a) 

7  •  (r  ,a)  =  oL  (r;w)/eo 
7  *  Hc  (r;a)  =  (r;a)/yQ 

Mote  that  the  time  dependence,  exp [j at},  has  been  sup¬ 
pressed.  The  total  electric  and  magnetic  field  intensities 
are  denoted  by  Et(r;u)  and  n~  ( r;w),  the  total  electric  and 
magnetic  current  densities  are  denoted  by  Jt(r,-w)  and 
5?T(r;a),  and  the  parameters  ?C(r;a),  m'(r;u),  e0,  Ua  ,  a,  and 
r  denote,  respectively,  totai  electric  and  magnetic  charge 
densities,  electric  permittivity,  magnetic  permeability, 
frequency,  and  observation  position. 


M  (r;a) 
a  (r;a) 


(2-1) 
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■r 


to 
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In  the  case  of  scattering  by  an  obstacle,  we  may  decom¬ 
pose  the  total  fields  and  sources  as 


p 

=t 

c. 

(r;oj) 

sine  ,-  . 

=  =■  (r;w)  + 

ES  (r;w) 

(2-2a) 

n 

=t 

a 

(r;w) 

=  (r;w)  + 

HS  (r ; lo) 

(2- 2b) 

3C 

(r;w) 

*  3s  (r;w)  +  J 

(?;w) 

(2-2c) 

MC 

(r;w) 

*  >iS  (r;a)  +  M 

(r;w) 

( 2— 2d) 

;« 

,  t 

(r;->) 

=  oS  (r;jj)  +  5 

(r;w) 

(2-2e) 

.  t 
a 

(r;w) 

*  aS  (r;w)  +  m 

(r;w) 

( 2— 2f ) 
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inc 

r  inc 

=  s  =  s 

where  E 
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the  incident  fields 

which  J  ,  M  , 

p  ° ,  and 

-  s 

m  g 

ive  rise  to. 

,  =  s  =.  s 

and  E  ,  n.  are  tne 

scattered  fields 

due  to  the  sources  J,  M,  p,  and  m  induced  on  the  scaCterer 

i'fcr  dielectric  scattering  these  scurces  are  interpreted  as 

effective  sources  that  replace  the  obstacle).  The  scattered 
-  s  -  s 

fields  S  ,  H  obey,  then,  the  vector  Helmholtz  equations 


7  x  7  x  E^  (r;w)  -  k“Es  (r;w)  =  -jwuQ  J  (r;u)  -  7  x  M  (rpx) 


(2-3a) 


and 


7  x  7  x  nS  (r;x)  -  k"  HS  (r;w)  =  -jueo  -1  ^7x3  (r;u) 

(2- 3b) 

where  k  is  the  wavenumber,  { s , u ,  )  . 

The  fields  which  satisfy  (2-3)  may  be  found  by  convolving 
the  impulse  response  cf  (2-3)  with  the  forcing  functions 


present.  The  impulse  response  is  obtained  d y  determining 
the  dyadic  Green's  function,  r(r,r';u)  (a  double  bar  will 
denote  a  dyadic),  which  satisfies 


x  7  x  r  (r,?;_i)  -  k“  T  (r,?;u)  *  5  (r.r*)  I 


(2-4) 


Here,  I  denotes  a  unit  dyadic,  r'  denotes  the  source  posi 
ticn,  and  5(r,r')  denotes  a  three-dimensional  Dirac  delta 
distribution.  The  solution  of  equation  (2-4)  is  given  by 


(r,?;o)  =  (I  77)  G  (r,?;u>) 


(2-5) 


where  G(r,r';w)  is  the  free-space  Green' 


s  function 


■jkI  r  -  r ' 


(2-6) 


—  S  —  —  q  « 

The  scattered  fields  E  (r;w)  and  h  (r;a)  may  new  be  ex¬ 
plicitly  represented  by 


Is  (r;w)  =  '  f  (r ,r ’  ;w)  •  ( —  J -»u Q  J  (?:u)  -  :t  M  (?;*)]  dr’  (--7a) 


(r:w)  =  ,  7  >1  "  x  J  (d;-)]  dr’  C-7b) 

zrr  f 


;here  V  denotes  the  volume  occuoieti  bv  th? 


We  substitute , 


next ,  equations  (2-7a,b)  into  equations 
(2-2a,b).  Using  the  vector  identities  G7 ' xM=7 ' x(GM) -7 ' GxM 
( 7 '  G )  •  7 1 xfl= - V ' • ( V '  Gx>! ) ,  and  the  relation  7G(r,r';w)= 

-V* G( r, r' ;u) ,  we  obtain  the  following  space-frequency  repr 
sentations  for  the  total  electric  and  magnetic  fields: 


at,-  N  =inc,-  .  ,  1 

&  (r;u)  =  c.  ( r ; oj) 


J 


[k"I+77]  •  j(r')G(r,r'  ;jj)dr' 


±  n-r.. 


r 


-—7  [k  1+77  •  1  i  n  xM(r  )G(r ,r  ;a)ds 

k2  !  >  ouC 

*  I  3 


M(r')  x  7*G(r ,r '  ;a»)dr '  (■  r  t  V' 


(2-8a) 


V' 


ana 


a’U;*)  »  +  7-7—  ;  [k“I+77]  •M(r')G(r,r’  ;w)drf 

J  ^  j 

0  v 


f 


+  -^7  rk”I+77 j  *  J  !  n  xJ(r’ )G(r,r '  ;^j)d£ 
.  z  out 

K  ;  s 


+  5  J(r')  x  7'G(r,?;^)dr'  >  r  i  V 


(2-Sb) 


The  integration  over  the  surface,  S,  denotes  integration 
over  the  surface  bounding  the  volume  V' . 

The  space-time  representations  of  the  electric  and  mag¬ 
netic  fields  may  be  obtained  by  inverse  Fourier  transformi 
the  frequency  dependence  found  in  expressions  (2-Sa)  and 
(2-Sb)  i 9 ' .  They  are  given  by 


The  representations  presented  are  general  expressions 


which  are  valid  for  an  arbitrary  scatterer  positioned  in  the 
previously  assumed  medium  for  all  r  such  that  r= r ' .  At  the 
offending  point  r=r ' ,  the  expressions  become  singular  and 
hence  must  be  evaluated  by  considering  the  limit  as  r  ap¬ 
proaches  r'  [9,12].  The  Cauchy  or  Kadamard  principal  value 
[13]  is  typically  used  for  the  description  of  these  integ¬ 
rals.  The  frequency-domain  representations  of  the  electric 
and  magnetic  fields  become  in  the  Cauchy  principal  value 
sense  (a  single  bar  through  the  integral  will  denote  a 
Cauchy  integral) 


St,- 


E‘"(r ;u)  -  2ElnC(r;u)  +  FP  [k2!+77]  •  J(r')G(r,r'  ;u)dr' 

^  "“"n  Jv 


2= 

--T  [’*  1+77] 

k“ 


r 


n  M(r,)G(r,r,;u)ds' 
5 

1 


-  -  M(r')  x  7  G(r jj) dr '  >  r  £  V' 

j 


(2-10a) 


and 


FP 


V' 


[k2 1+77 ] 


M(r,)G(r,r,;w)dr‘ 


[k"I-: 


x  J(r')Gl.r,r'  ;w)ds' 
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Note  that  an  interchange  of  primed  to  unprimed  coordinates 
has  been  made.  Similar  factors  of  two  appear  in  the  time- 
domain  representations,  and  F?  denotes  'the  finite  part  of'. 

These  results  may  be  specialized  to  describe  thin,  per¬ 
fect  conducting  surfaces  [12].  On  such  a  surface,  the  ap¬ 
propriate  boundary  conditions  [10]  are  that  the  tangential 
total  electric  field  is  zero,  i.e.,  nxE  =0  (n  defined  to  be 
the  outward  normal  unit  vector  on  S),  and  that  the  tangen¬ 
tial  total  magnetic  field  is  equal  to  an  equivalent  surface 
current  source,  J  ,  i.e.,  nxHt=Js .  With  these  boundary  con¬ 
ditions,  we  may  immediately  write  the  space-f requency  repre¬ 
sentations  for  the  electric  and  magnetic  fields  on  the  sur¬ 
face  S  as 


;inc  -  . 

-n  x  a  (,r;ui) 


nvFP  •  2  ~  ~ 

ISTi  [kI*  ”!•*.<?’>•>>  G  <*•*’>“>  (Ml.) 

o  '  s 


and 


j  (r;j)  =  2n  x  H"nc(r;jj)  -  2nx-  J  (?;u)  x  7G(r  ,r 1  ;^)  ds '  (2-llfc) 
s  Js  3 


3imi*ar^.y,  tne  space-tame  representations  are  given  by 


:  ,-iac  , -  n  _ ,  - -l  :  _ 

:0  n  X  5t  L  vr;C'  *  ir  X  ^  2  1  ' j  x 

'  S  C~':Z~ 


( 2-12a) 


T  =  C-R/c 
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and 


J  (r;t)  =  2  n  x  Hlnc(r;t)  +  x  ~ 
s  2lT  J  S 

+  4-1  (r';T)  x  (r~r'}]  ds’ 

3T  3  cR2 

Note  that  the  equivalent  magnetic 
not  appear  in  these  expressions, 
electric  field  by 


[Js(r';T)  x 


(r-r') 


T=t-R/ c 

surface  current, 

R  is  related  to 
s 


(2-12b) 

M  ,  does 

5 

the  total 


M 


-n  x  £ 


(2-13) 


which  vanishes  for  perfect  conducting  surfaces.  Note,  also, 
that  for  good  conductors  the  effective  current  source  J  may 
be  replaced  bv  cl1  (o  denotes  the  conductivity  of  the  obsta¬ 
cle),  and  therefore  terms  involving  nxJ  also  tend  to  zero. 

As  a  final  remark,  we  note  that  the  term  2nxHinc  appear¬ 
ing  in  expression  (2-12b)  is  commonly  known  as  the  physical 
optics  approximation  for  the  current  density  Js  .  This  ap¬ 
proximation  is  useful  for  testing  the  validity  of  results 
obtained  from  expressions  (2-12a,b)  when  no  results  for  com¬ 
parison  exist. 


Chapter  III 

FUNDAMENTAL  CONCEPTS  OF  THE  SINGULARITY 
EXPANSION  METHOD 


3 . 1  INTRODUCTION 

The  motivation  for  the  singularity  expansion  method  (SEM) 
is  essentially  based  on  experimental  observations  which  have 
established  that  the  transient  surface  currents  generated  on 
structures  (scatterers)  by  arbitrary  excitation  are  primari¬ 
ly  in  the  form  of  damped  sinusoids;  the  particular  shape  be¬ 
ing  dependent  on  the  form  of  excitation  and  the  specific 
geometry  of  the  structure  under  consideration.  3y  assuming 
that  a  scatterer  can  be  uniquely  specified  mathematically  by 
an  associated  modal  and  pole  structure,  and  that  the  form  of 
the  excitation  is  known,  the  SEM  enables  one  to  determine 
the  surface  currents  directly  from  an  appropriate  expansion 
of  these  parameters.  Specifically,  the  expansion  was  found 
to  require  knowledge  of  four  parameters  [I]:  the  natural 
frequencies  and  corresponding  natural  modes,  the  structure 
of  the  incident  wave,  and  scalar  coefficients  that  couple 
the  natural  resonances  to  the  incident  wave  (coupling  coef¬ 
ficients).  Since  the  form  of  the  excitation  is  assumed  tc 


ce  ,<n own,  tne  nature,  i recuer.cies ,  r.atura.  mooes,  ar.c 


r.g  coefficients  need  to  be  determined  i; 
h  an  SEM  representation  of  the  problem. 


oroer  to  estao- 


Ma tnemati c ally,  the  expansion  for  the  space- frequency 
surface  currents  induced  by  delta  function  excitation  on  fi¬ 
nite,  perfect  conducting  objects  in  free  space  is  given  by 


U(r;s)  =  >  n  v  (r)  (s-s  )  +  W(r;s) 

Ct  Ct 

a 


In  the  time-domain,  this  recresentation  becomes 


U(r;t)  -  £  na  v  (r)  2  +  W(r;t) 

cl 


(3-la) 


(3-lb) 


In  these  equations,  s  is  a  complex  variable  which  is  related 
to  the  frequency,  w,  by  Im{sl=u,  U( r ; s ) ,  U(r;t)  denote  the 
space-frequency  and  space-time  surface  currents,  n.  denotes 
the  coupling  coefficient  associated  with  the  pole  s. ,  v  (r), 
va(r)  denote  the  natural  mode  vectors  associated  with  s  , 
W(r,*s)  denotes  an  entire  function  and  W ( r ;  t )  denotes  the 
corresponding  time-transformed  function,  m^  denotes  the  mul¬ 
tiplicity  of  the  pole  sa ,  and  the  summations  are  over  all 
poles.  In  Section  3.2.1,  we  consider  space-frequency  tech¬ 
niques  for  obtaining  the  natural  frequencies,  and  natural 
modes.  In  Section  3.2.2,  we  present  space- frequency  techni¬ 
ques  for  obtaining  the  coupling  coefficients,  and  briefly 


Section  3.3.1,  we  develop  the 
Icrdaro-Davis  method  for  obtair.incr  the  natural  resoonses. 
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Available  techniques  fcr  analyzing  the  stability  of  various 
finite  difference  approximation  schemes  are  also  discussed 
in  that  section.  And  in  Section  3.3.2(  we  present  transient 
matrix  methods  for  determining  the  coupling  coefficients. 

3.2  SPACE - FREQUENCY  TECHNIQUES 
3.2.1  Natural  Erecuencies  and  Modes 

An  arbitrary  Fredholm  integral  equation  of  the  first  kind 
(e.g.,  expression  (2-lla))  may  be  cast  in  the  general  form 

r  ^ 

j  r(r,r';s)*U(r,;s)dr'  =I(r;s)  (3-2) 

R3 

where  r  (  r ,  r ' ; s )  denotes  a  dyadic  kernel,  U ( r ' ; s )  denotes  the 
desired  unknown,  and  I(r,-s)  denotes  an  arbitrary  forcing 
function. 

For  simplicity,  we  will  write  these  integral  equations 
using  the  inner  product  notation  [ 1 ] 

<r(r,r';s);  U(r';s)>  =  i(r;s)  (3-2a) 

where  the  appropriate  operation  between  the  kernel  and  unk¬ 
nown  will  be  given  above  the  comma  separating  these  parame¬ 
ters,  and  the  integration  is  with  respecc  to  the  common  spa¬ 
tial  variable. 


A  natural  mode,  v, (r),  satisfies  equation  (3-2a)  in  the 
absence  of  a  forcing  function.  We  may  write 


<f(r  ,r  '  ;s^) ;  .^(r')>  =  0 


(3-3) 


where  s~  denotes  the  corresponding  complex  natural  frequen- 


The  parameters  s^  and  v2(r)  may  be  found  by  discretizing 
equation  (3-3)  using  a  method  of  moments  (14]  formalism.  We 
obtain 


(f  (s  )}  •  (V  )  -  (6  ) 

^  n,m  a  ;  n  a  n 


(3-4) 


where  n,  m  are  positive  integers,  (F  _  (  =  „))  denotes  an  r.  bv 
m  matrix,  (v^)^  denotes  the  unknown  mode  vector  of  length  n, 
and  (0n)  is  a  zero  vector  of  length  n.  The  magnitude  of 
both  n  and  m  is  dependent  on  hew  refined  the  discretization 


Equation  (3-4)  represents  a  homogeneous  sysrem  of  equa¬ 
tions.  Such  a  system  has  a  solution  if  and  only  if  the  ma¬ 
trix  ( f  n  (s2))  is  singular.  Hence,  the  natural  frequen¬ 


cies,  sa,  may  be  found  by  solving 


det  (F  (s  )]  =  0 
n,3i  3 


(3-5) 


The  natural  modes  may  now  be  found  from  equation  (2-4)  using 
the  results  of  equation  (3-5). 

Equation  (3-5)  is,  in  general,  extremely  complicated  to 
solve.  Numerical  solution  techniques  typically  use  either  a 
function  iteration  root  searching  technique  (Section  4.2.1) 
or  a  contour  integration  [15]  (Section  4.2.2)  method.  The 
use  of  contour  integration  allows  one  to  locate  desired 
roots  by  partitioning  the  complex  plane. 


3.2.2  Coupling  Coefficients  and  Entire  Functions 

The  following  derivation  for  obtaining  the  SEM  coupling 
coefficients  patterns  a  development  due  to  Baum  [5]. 

Associated  with  the  coupling  coefficient,  6  .,  is  a  cou¬ 
pling  vector,  ^(r).  The  coupling  vector  is  defined  to  be 
the  conjugate  adjoint  of  the  natural  mode,  va(r),  and  hence 
satisfies 

<u^ (r ' )  ;  *(?,?' ;sa»  =  0  .  (3-6) 


3y  applying  the  method  of  moments,  we  have 


(-  >  *  0*  (S  )]  =  (0) 

n  a  v  rum  a  n 


(3-7) 


The  kernel  is  now  expanded  in  a  Tavlor  series  about 


(3-3a) 


r(r,r ' ;s)  =  /  (s-s  ) "  f,  (r,r’) 

i=0  a 


=  fr^r 

"  e*  1  «?  =  c 


( 3— 3b ) 


Che  forcing  function  is  similarly  expanded  as 


I(r;s)  =  l  (s-s  )  1  (r) 

i=0  01 


(3-9) 


’  TT^T  !<;;s>|- 

as  ! S=S 


Assuming  only  a  first  order  pole,  we  may  write  the  res¬ 


ponse  from  equation  (3-la)  as 


Z (r ;s)  = 


V"ct(r)(S-S-x) 


13-10) 


where  U'(r;s)  denotes  some  analytic  function  about  s=s-,. 

3y  substituting  (3-9),  (3-10)  and  (3-11)  into  the  basic 

equation  (3-2a)  and  matching  powers  of  (s-s,)  ,  we  obtain 


<:  . (r,r');  n  v  (r)>  =  0 

U ,  ji  j.  j. 


(3-lla) 


<?.  (r , r ' ) ;  * 


;^(r)>  -  <f  (r.r');  U'(r':s)>  =  I  (: 


(3-11 b) 


Operating  on  (3-1 2b)  from  the  left  by  u  (r)  yields 


I 
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<:,(r):  "l,a(r’r'):  V»(r)>  =  <:,(r);  !o,a(r)> 


(3-12) 


sir.ce 


<u  (r) ;  ?  (r , r ' ) ;  U'(r';s)>  = 

Qt  U  ,  Cl 


f  :  r  _ 

dr'  I  dr  »  (r) 

J  J  * 

_3  L  _3 

R  R 


(r,r')|  •  U' (r* ;s)  =  0 


(3-13) 


by  equation  (3-6).  Therefore, 


<£,(*>;  ~l0  a(r)> 

- -  (3-1.4) 

;:..(r);  ",  (r ,  r ' ) ;  ■'  (?)> 


is  the  expression  for  the  coupling  coefficient  ai  s=:s2- 

The  coupling  coefficients  relate  the  incident  waveform  to 
che  modal  structure  of  an  object.  They  indicate  which  modes 
are  excited  and  the  extent  to  which  they  are  excited.  3aum 
[5]  has  discussed  two  different,  but  ultimately  equivalent, 
types  of  these  coupling  coefficients  in  order  to  treat  two 
different  philosophical  interpretations  as  to  how  modes  are 
activated.  In  one  interpretation,  all  modes  are  excited 
simultaneously  across  an  object  no  matter  where  on  the  ob¬ 


ject  tr.e  excitation  oricrmaced. 
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cion,  modes  in  various  regions  cannot  be  excited  until  the 
incident  wave  has  reached  those  regions.  We  will  not  pursue 
these  types  further  here. 

The  entire  function,  W,  associated  with  the  pole,  s  ,  is 
necessary  for  equation  (3-la)  to  be  mathematically  valid. 

Its  form  and  use  are  not  well  understood,  however.  Typical¬ 
ly,  the  entire  function  is  omitted  by  the  empiric  justifica¬ 
tion  of  obtaining  current  distributions  directly  from  a  set 
of  poles  which  are  in  good  agreement  with  the  distributions 
obtained  by  standard  methods  [5,16].  The  physical  signifi¬ 
cance  of  the  inclusion  or  omission  of  the  entire  function 
requires  further  consideration. 

3.3  SPACE-TIME  TECHNIQUES 

3.3.1  Natural  F recuencies ,  Natural  Modes ,  and  Stability 
Considerations 

In  the  time-domain,  electromagnetic  integral  equations  of 
the  first  kind  may  be  written  in  general  form  as 

r  ( 

|  I  “(r.r'jt-t')  •  UCr'jt'jdt’dr'  =  X(r ; t)  (3-13) 

J  > 

3  1 

R  R 

where  F( r, r' ; t-t' )  denotes  a  retarded  dyadic  Green's  func¬ 
tion  [11],  U(r';t)  denotes  the  desired  unknown,  and  I(r,-t) 
denotes  an  arbitrary  forcing  function. 


For  illustrative  purposes,  we  will  restrict  the  discus¬ 
sion  in  this  section  to  thin,  perfect  conducting  surfaces 
for  which  integral  expression  ( 2-12a)  is  appropriate.  The 
discussion  will  also  be  limited  to  rectangular  ( x , y , 2 )  coor¬ 
dinate  systems.  A  similar  development  applies  to  other  ex¬ 
pressions  which  may  be  cast  in  the  form  of  equation  (3-15), 
and  other  coordinate  systems. 

Since  the  spatial  differential  operators  appearing  in  ex¬ 
pression  (2-12a)  are  with  respect  to  the  unprimed  coordi¬ 
nates,  the  following  variation  of  this  expression  is  valid: 


_L  r^nc 

r  C 


C) 


n 


J  (r * ;  c-R/c) 


-7T 


dr ' 

(3-16) 


Here,  n  is  the  outward  normal  on  some  arbitrary  surface  S. 
The  integral  over  this  surface  is  commonly  known  as  the  mag¬ 
netic  vector  potential.  By  letting  A(r;t)  denote  this  po¬ 
tential,  we  may  write  (3-15)  as 


O 


The  current  density,  J  ( r ; t-R/c ) ,  aocearing  in  (3-15)  is 
typically  the  unknown  which  is  desired.  However,  for  r.ota- 
tional  purposes,  and  stability  analysis,  expression  (3-15a) 


grass . 


4 


r 

V 
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In  general,  when  a  desired  unknown  appears  buried  within 
the  integrand  of  an  integral  equation  it  is  not  possible  to 
determine  it  analytically.  To  obtain  a  numerical  solution, 
one  generally  begins  by  expanding  this  unknown  in  some  sui¬ 
table  set  of  basis  functions.  If  the  function  to  be  expand¬ 
ed  is  at  least  piecewise  continuous  over  the  region  of  in¬ 
terest,  a  suitable  basis  set  would  be  a  pulse  expansion  for 
the  spatial  variables.  If  the  function  is  also  reasonably 
well  behaved  through  time,  the  temporal  dependence  may  also 
be  expanded  as  pulses.  The  function  Js(r;t-R/c)  generally 
satisfies  these  requirements,  and  hence  a  pulse  expansion  in 
both  space  and  time  is  appropriate.  It  should  be  noted, 
however,  that  this  approximation  can  become  quite  peer  at 
surface  edges  due  to  the  singular  behavior  of  the  current 
component  parallel  to  the  edge.  Special  care  is  required 
for  such  structures  (Chapter  5). 

The  expansion  of  the  current  density  may  be  written  as 


J„ (r ; t-R/c) 


N 

~  J.  ?,  (t-oat-R/c) S . (r ) 
-  1?  at  1 

1=1  p=-=° 


(3-17) 


where 

(  1,  for  t  in  the  time  interval  centered  at  p  It-tR/c 
?  ^(r-pic-R/c)  =  J 

i_  0,  elsewhere 


a 


r) 


■  1,  for  r  in  the  space  segment  centered  at  ilr 

t 

'l 

i  0,  elsewhere  . 


l~s 


Here,  J.  denotes  the  current  amoiitude  coef f icients ;  i 

i,? 

denotes  a  general  spatial  index,  i.e,  i  may  represent  one, 
two  or  three  integer  variables  depending  on  the  geometry  of 
the  problem;  N  denotes  a  general  upper  bound  for  the  surarca 
tions  corresponding  to  each  of  integer  variables  which  i 
represents;  and  Ar  denotes  a  general  spatial  sampling  dis¬ 
tance,  i.e.,  A r=  (Ax, Ay, At). 

Expansion  (3-17)  enables  one  to  write  the  vector  poten¬ 
tial  appearing  in  expression  (3-16a)  explicitly  as 


A(mA,nA,kA;pAt) 


N  N  N, 
m  n  k  =° 

:  I  I  l  Ja.  n.  k.  P. 

m'-l  n'-l  k'-l  ’  ’  ’? 


,  n-n*  ,  k-k' 


where 


(.1+1/2)  A  (3+1/  2)  A  ( -2+1  / 2) A 


?  ,  _  (n-(u“+v‘'+w")  '  -/c) 


'  !  i  i 

(ci-1/2)  A  (3-l/2)A  (5-1/ 2) h  ^  ; 


dudvdw  , 


m,  n,  and  k  being  positive  integers  which  are  bounded  by  N 
Nn,  and  Nk  respectively,  and  p  is  an  unbounded  integer  (by 
causality,  p  may  be  restricted  to  positive  integers).  Not 
that  in  this  expansion  we  have  tacitly  assumed  that  the  sp 
tial  sampling  distance  is  uniformly  equal  to  some  constant 


sc  that  the  continuous  variables,  (x,y,z),  correscor.d  t 


the  discrete  variables,  ( mA ,  r.A ,  kA )  .  This  is  typically,  but 
not  necessarily,  done. 

The  current  coefficients  appearing  in  (3-18)  are  the  de¬ 
sired  parameters.  They  may  be  extracted  by  approximating 
the  continuous  differential  operators  appearing  in  (3- 15a) 
by  central  finite  difference  operators.  A  thorough  discus¬ 
sion  of  finite  difference  approximation  (FDA)  techniques  may 
be  found  in  Ames  [ 17 ] . 

In  passing,  it  is  worthwhile  to  note  that  it  is  possible 
to  establish  an  analytic  equivalence  between  the  finite  dif¬ 
ference  formulation  of  a  time-domain  problem  and  the  basis 
set  formulation  of  the  equivalent  frequency-domain  problem 
by  using  inverse  transform  techniques.  This  equivalence  is 
satisfying  since  it  establishes  that  finite  difference  tech¬ 
niques  are  not  simply  convenient  mathematical  tools  for  the 
solution  of  time-domain  problems,  but  are  appropriate,  phy¬ 
sically  meaningful,  methods  of  solution. 

By  using  finite  differences,  the  time  derivative  of  the 
vector  potential  may  be  written  as 


“T A(r;t)  =  — - — -  [A(r ;  (?+i)Ac)  +  A(r :  (p-l)it) 
c"  3 t“  (cit)" 

-  2  A(r ;pit) I  +  Or(it)2' 


(2-19) 


where  0((At)7-)  denotes  the  order  of  the  truncation  error  in¬ 
troduced  in  the  FDA.  The  soatial  ccerators  may  be  similarlv 

-  -  a  * 


differenced  (Section  5.3,  and  Sections  5.3,  6.4] 


Expression  ( 3-lSa)  nay  now  be  written  as  (p=l ,2 , . . . ) 


n  *  A(r ;  (p+l)iit)  =  (c±c)“  [n  *  7(7.  A(r;pAt) )  ] 


n  *  A(r;  (p-l)<lt)  +  2n  *  A(r;pAt)  + 


0(cit)“n  x  4 J—  Elnc(r:pAt)  +  (0(£t )"] 


(3-20) 


This  formulation  establishes  an  explicit  or  time-marching 


finite  time-differenced  scheme  for  the  vector  potential.  An 
explicit  scheme  allows  one  to  find  future  values  in  terms  of 
previous  results  without  the  need  for  a  matrix  inversion. 
Note  that  the  values  of  the  vector  potential  at  two  previous 
times  are  recuired. 


ferer.ce  operators,  and  by  using  expansion  (3-18)  to  repre¬ 


sent  the  vector  potential,  the  following  general  representa¬ 
tion  for  the  current  density  coefficients  is  obtained  (note 
that  the  notation,  TD+1  ,  used  to  represent  all  discrete 
functions  is  interpreted  as  T(mA,  nA,  kA,  (p+1 )  At )  )  : 


where  3D.  denotes  coefficient  matrices  corresponding  to 
different  times  {3_1  is  a  diagonal  matrix  corresponding  to 
J'p+1);  ?p+1  denotes  the  forcing  function  at  the  (p+l)-th 
time  step;  and  NT  denotes  an  integer  which  is  one  fewer  than 
the  number  of  time  steps  required  for  a  wave  to  propagate 
across  the  maximum  distance  of  the  structure;  in  other 
words,  if,  for  example,  six  time  steps  are  required  for  a 
wave  to  travel  this  maximum  distance,  NT  would  be  five  since 
the  summation  begins  at  zero.  The  prime,  indicates  a 

vector  of  the  currant  density  coefficients  of  every  spatial 
point  of  interest  on  the  structure.  And  as  a  final  remark, 
we  note  that  the  rank  of  the  3  matrices  is  dependent  on  the 
particular  geometry  of  the  problem  being  studied.  Tor  con- 


« 
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venience,  we  define  the  rank  of  these  matrices  no  be  seme 
integer,  N. 

To  obtain  the  natural  frequencies  and  natural  modes,  we 
are  interested  in  equation  (3-21)  in  the  absence  of  a  forc¬ 


ing  function,  i.e.. 


C  ,  J'  , 
?  P-P 


(3-22) 


where  C  ,=(3  ,  )  B  , .  The  solution  of  this  difference  ecua- 
P  “1  P 

tion  may  be  obtained  by  z  transform  techniques.  For  simple 
poles,  the  solution  is  given  by 


ti  p^l  - 

J '  =  z  V 

p+i  i  a 


(3-23) 


where  z^=exp  { s^  it  1  denotes  the  transient  representation  of 
the  pole,  s  ,  and  v3  denotes  a  vector  spatially  describing 
the  natural  mode.  For  poles  of  multiplicity  m„_, ,  *  1 ,  the 

solution  is  given  by 


(J*  )  -  (p+1)  z?  -j 

p-^-1  :t  a, a. 


(3-24) 


where  v ,  denotes  the  natural  mode  vector  corresponding  to 
a  pole  of  multiplicity  m2.  Note  that  entire  functions  do 
not  appear  in  this  development;  a  pole  structure  only  is  the 
basis  for  this  method.  Pole  clusters  mav  athsmst  to  model 


an  en: 


'unction  however,  and  therefore  entire  functions 
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nay  still  be  significant  although  they  are  not  explicitly 
represented  in  the  formalism. 

By  substituting  equation  (3-23)  (assuming  first  order 
poles)  into  (3-22),  we  obtain 

[—  N^.  -j 

I  -  l  C  ,  z'(p,+1)  v  =  0  .  (3-25) 

p'-O  ?  01  j  a 


This  is  a  homogeneous  system  of  equations, 
found  from 


det 


N  1 

I  -  l  C  ,  z-<”+1)!  -  0 

L  p'-o  p  ’  J 


The  poles  may  be 


(3-26) 


The  modes  may  now  be  found  from  (3-25). 

There  is  an  alternative  to  this  z  transform  solution 
technique.  Any  finite  difference  scheme  in  the  form  of 
equation  (3-22)  may  be  condensed  into  an  equivalent  two- lev¬ 
el  matrix  form  [17]  by  introducing  a  state  vector,  for 
the  p-th  time  step,  such  that  (T  denotes  transpose) 


-T 

K 

D 


L_ 


(3-27) 


transition  matrix,  f,  such 


a  state 


that 
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N 


M 


i  = 


co  ci  .  % 

I  0  .  0 

0  1. 

•  •  •  • 

•  •  •  •  • 

o  .  o  i  o 


(3-28) 


and  forming 


K  . .  =  ?K 
P+1  P 


(3-29) 


A  discussion  of  error  propagation,  or  a  stability  analy¬ 
sis  of  finite  difference  schemes  is  appropriate  at  this 
time.  A  discrete  finite  difference  representation  of  a  con¬ 
tinuous  problem  may  yield  an  unstable  (unbounded)  solution 
when  certain  relationships  between  the  sampling  distances 
used  for  different  variables  are  not  satisfied.  For  hyper¬ 
bolic  equations  (wave  equations,  e . g . ,  equation  (3-16))  the 
relation  between  the  time  sampling  (At)  and  spatial  sampling 
(A,  assuming  a  uniform  sampling  distance  in  all  directions) 
distances  are  of  interest.  It  has  been  shewn  by  Courant, 
Friedrichs  and  lewy  (CFL)  [17]  that  the  time  sampling  dis¬ 
tance  for  these  equations  can  be  at  most  equal  to  the  spa¬ 
tial  distance,  i . e . ,  At=A .  This  is  the  most  lax  restriction 
possible;  it  can  tighten  consiaerablv  demanding  on  how  the 


a 
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discretization  is  implemented.  Two  methods  are  available  to 
analyze  the  stability  of  finite  difference  schemes  for  li¬ 
near  equations.  We  consider  these  now. 

The  stare  transition  matrix  appearing  in  equation  (3-2S) 
contains  ail  the  information  of  the  finite  difference  ap¬ 
proximation  (including  boundary  conditions).  A  stability 
analysis  of  the  difference  scheme  may  be  done  by  examining 
the  magnitude  of  the  eigenvalues  of  this  matrix.  If  all  the 
eigenvalues  are  less  or  equal  to  one  in  magnitude,  errors 
will  not  grow  through  time  and  hence  the  solution  will  be 
bounded.  This  technique  is  known  as  matrix  stability  analy¬ 
sis.  The  matrix  stability  method  is  useful  for  testing  if  a 
known  CFL  condition  yields  a  stable  solution.  It  does  not 
predict,  in  general,  the  specific  numerical  value  required 
for  stability.  An  alternate  method  may  be  used  to  deter¬ 
mine,  or  at  least  approximate,  this  value. 

A  simple  method  known  as  Fourier  stability  analysis  may 
be  used  to  determine  the  stability  criterion  for  an  uncom¬ 
pressed  difference  scheme  (e.g.,  equation  (3-20)  or  (3-22) 
instead  of  equation  (3-29)).  The  method  analyzes  only  the 
specific  difference  equation  and  hence  ignores  the  influence 
of  boundary  conditions.  Since  boundary  conditions  can  in¬ 
fluence  the  stability  of  a  scheme,  the  Fourier  method  is  not 
considered  as  thorough  as  the  matrix  method.  However,  since 


a  specific  number,  whether  exact  cr  approximate,  for  the  CFL 
condition  is  readily  produced,  this  method  provides  useful  a 
priori  information  about  a  particular  difference  formula¬ 
tion.  The  matrix  method  may  always  be  used  to  confirm  the 
stability  criterion  given. 

In  brief,  the  Fourier  method  examines  the  propagating  ef¬ 
fect  of  a  single  row  of  errors  along  some  arbitrary  line  of 
the  FDA.  This  is  accomplished  by  determining  an  exponential 
solution  for  the  difference  scheme  from  discrete  separation 
of  variable  techniques.  For  a  stable  solution,  restrictions 
on  the  exponential  solution  must  be  enforced.  A  one-dimen¬ 
sional  example  may  be  found  in  Section  5.3.  Two-dimensional 
examples  may  be  found  in  Sections  5.3,  5.4. 

Stability  alone  does  not  imply  convergence  of  the  FDA  to 
the  true  solution.  For  a  thorough  discussion  on  matrix  and 
Fourier  stability  methods  and  convergence  requirements  one 
should  refer  to  Ames  [17]. 

The  stability  of  physical  problems  is  mathematically  de¬ 
scribed  by  the  location  of  poles  in  the  complex  plane.  The 
stability  of  the  finite  difference  representation  of  elec¬ 
tromagnetic  expressions  is  dependent  on  the  magnitude  of  the 
eigenvalues  of  the  state  transition  matrix.  Hence,  we  anti¬ 
cipate  some  relation  to  exist  between  these  eigenvalues  and 
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The  relation  follows  simply  by  considering  the  solution 
of  the  difference  equation,  equation  (3-24),  applied  to  the 
state  transition  formulation.  For  first  order  poles,  we  may 
write 


;?+1  K  =  i  z?  K 
ct  ct  ct  ct 


(3-30) 


2  K  »  0  K  . 
a  a  a 


(3-31) 


This  represents  an  algebraic  eigenvalue  problem  for  the  ei¬ 
genvalue,  and  the  eigenvector,  K  .  It  can  be  shown  that 


the  natural  mode,  v  ,  and  X  are  related  to  one  another  by 

i  J. 


—  \T  M 

-T  i  T  -T  T  -  T 
K  =  j  z  v  ,  z  v 
a  a  a  a 


(3-32) 


mhe  poles,  s^,  may  be  found  by  solving 


s  =  ln(z  )/At 
a  a 


(3-33) 


3.3.2  Coupling  Coefficients 

A  method  has  been  presented  which  determines  the  natural 
frecuencies  and  natural  modes.  To  complete  the  SEM  form  of 


solution  we  need  to  determine  the  coefficients  that  coup!; 


me  natural  frequencies  ana  .T.oaes  no  me  incident  forcing 
function.  Two  different  formulations  of  these  coupling 
coefficients  are  possible  in  the  time-domain.  Cordaro  [13] 
has  suggested  a  method  to  obtain  an  exact  representation  of 
the  coefficients  when  a  complete  set  of  distinct  eigenvalues 
(first  order  poles)  is  known.  This  is  accomplished  by  using 
eigenvector  decomposition  techniques.  The  basic  method  may 
be  extended  to  obtain  an  approximation  to  the  coefficients 
when  only  a  partial  set  of  distinct  eigenparaneters  is 
known.  An  alternate  formulation  for  a  partial  set  of  poles 
is  a  time-domain  analog  of  the  f requency-dcmain  technique 
previously  presented  (Section  3.2.2).  We  will  initially 
consider  Cordaro 's  method. 

Define  a  state  vector,  U3  ,  to  represent  a  normal  incident 
forcing  function  at  the  p-th  time  step  as 

"T  =  [F  ,  0 . OP  .  (3-34) 

?  P 

Here,  F  and  C  are  2xN  row  vectors. 

p 

Note  that  for  a  delta  or  impulse  excitation  only  U0  is 
nonzero.  The  state  current  distribution  may  new  be  written 
explicitly  as  (a  state  representation  of  the  forcing  func¬ 
tion  has  been  added  to  equation  (3-29)) 


t 
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c  =  r"1  U0  . 


(3-40) 


Combining  these  results  yields 


T  =  MT  Ap_1  C 


(3-41) 


or  equivalently 


—  r  0—1  — 

T  =  )  z‘  n  v 
p  u  at  ct  a 
a 


(3-42) 


where  the  pole  sa  is  related  to  z3  by  equation  (3-33),  na 
denotes  the  coupling  coefficients,  and  denotes  the  cor¬ 
responding  natural  modes.  Equation  (3-42)  is  the  desired 
SEM  representation. 

When  only  a  partial  set  of  eigenvalues  is  known,  the  de¬ 
composition  of  $  as  given  by  equation  (3-36)  is  not  possible 
exactly  since  the  inverse  which  appears  only  exists  in  a 
generalized  or  pseudo  inverse  sense.  Therefore,  only  a 
least  squares  approximation  to  the  coupling  coefficients  is 
possible  in  this  case. 

This  complication  may  be  avoided  by  developing  a  time-do¬ 
main  formulation  for  the  coupling  coefficients  analogous  to 
the  frequency-  domain  method  (Section  3.2.2).  We  begin  by 
replacing  r ( r , r ' ; s )  by  a  matrix  function  T(z)  defined  by 


T(s)  =  2  " 


N„+i  ‘J  OT„-i) 


I 


3y  following  similar  power  series  expansions,  we  “her.  defi 
Tx  (za)  in  analogy  with  ?1>a  (r,  r' )  to  be 


V  (nt;l) 


T,  ,(*,)  =  j  (N„+l)z  I  -  (N  -i)C  z 

l'a  *  I  "  i=0  T  1 


(NT-i-D 


We  define,  next,  the  coupling  vector,  vta(r),  to  be  the  fir 
block  of  N  elements  of  the  left  eigenvector1  corresponding 
to  the  ex-th  eigenvalue  of  the  state  transition  matrix. 

3y  replacing  the  frequency-domain  inner  product  opera¬ 
tions  by  matrix  multiplications,  we  may  write  the  coupling 
coefficient  at  z=za=exp { s^ At }  as 

_T  - 

-■  (z  )  =  - - -  (3-3 

-  *  a1*!.  (2  )  V  J 

ct s  1  ,a  3  a' 


where  ^z^)  is  the  forcing  function  vector  evaluated  at  z.x 
It  should  be  noted  that  when  each  of  the  sub-matrices  o 
the  state  transition  matrix  are  symmetric,  the  first  N  row 
of  the  left  and  right  eigenvectors  are  identical  to  a  nor¬ 
malization  factor.  This  is  not  true  for  the  remaining  por 
tion  of  these  vectors,  however,  since  it  can  be  shown  that 


*  Let  p  oe  the  right  eigenvector  or  tne  transpose  of  some 
matrix^ A  corresponding  to  t h e ^ e i_g e ny a  1  u e  X.  Then  p  sat: 
fies  A*p=\p.  Now  consider  (A*p)~  =p~  A=p~  \  .  In  this  case 
p  is  known  as  the  left  eigenvector  of  the  matrix  A. 
Hence,  p  is  either  the  right  eigenvector  associated  with 
the  matrix  A1  or  the  left  eigenvector  associated  with  th 
matrix  A  corresponding  to  the  eigenvalue 


33 


the  left  eigenvectors  have  a  much  mere  complicated  structure 
than  the  right  eigenvectors. 

The  matrix  decomposition  formulation  is  recommended  when 
a  full  set  of  distinct  eigenparameters  is  known;  the  left- 
right  eigenvector  formulation  is  recommended  when  a  partial 
set  is  known. 

In  conclusion,  we  note  that  although  the  TD-SEM  formula¬ 
tion  for  obtaining  the  natural  frequencies  and  natural  modes 
is  relatively  straightforward,  a  fundamental  complication 
does  underlie  the  method.  Since  the  size  of  the  transition 
matrix  is  highly  dependent  on  the  geometry  and  the  level  of 
discretization  of  a  particular  problem,  it  is  possible,  even 
for  simp.e  geometries,  to  generate  a  transition  matrix  which 
surpasses  the  high  speed  storage  capabilities  of  the  largest 
computers.  A  variety  of  techniques  which  attempt  to  handle 
this  complication  by  taking  advantage  of  the  form  of  this 


matrix  are  oresented  in  ChaDter  4. 


Chapter  IV 

EIGENSOLUTION  METHODS  FOR  THE  TRANSITION  MATRIX 


4.1  INTRODUCTION 

The  TD-SEM  model  is  a  straightforward  and  efficient  meth 
od  for  determining  the  SEM  parameters  for  simple  geometries 
discretized  with  relatively  few  unknowns.  This  is  accom¬ 
plished  by  transforming  the  pole  searching  problem  into  an 
algebraic  eigenvalue  problem  (Section  3.3.1).  As  the  numbe: 
of  unknowns  increase,  however,  the  matrix  which  TD-SEM  gen¬ 
erates  becomes  unmanagably  large,  thereby  making  the  search 
for  eigenvalues  difficult  and  complicated. 

The  matrix  ?,  whose  eigensoiution  is  sought,  is  given  by 
equation  (3-28).  Some  comments  are  in  order  about  the  form 
and  properties  of  this  matrix. 

f  is  known  as  a  sparse  matrix  since  it  contains  a  large 
number  of  zero  elements.  It  is  in  block  upper  Hessenburg, 
or  more  specifically,  block  Frobenius  form  [19].  A  matrix 
in  Frobenius  form  possesses  no  symmetry  properties,  and 
therefore,  §  unfortunately  falls  into  the  class  of  unsymme- 
tric  real  matrices,  or  general  real  matrices.  This  is  in¬ 
deed  a  complication  since  the  field  of  eigensolution  method, 
is  both  narrowed  and  complicated  for  unsvmme trie  matrices 
due  to  the  possibility  of  obtaining  complex  eigenvalues  and 


genera-ized  eigenvectors. 
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A  matrix  in  Frobenius  form  possesses  the  property  that  it 
is  its  own  companion  matrix.  In  other  words,  the  problem  of 
determining  the  eigenvalues,  X,  of  4  may  be  dene  in  either 
of  two  possible  forms.  First,  we  may  consider  the  full  ma¬ 
trix  and  solve 

det  [4  -  XI]  =  0  (4-1)  , 

where  I  is  the  identity  matrix;  alternatively,  we  may  solve 

dec  (\NT  I  -  A^C  -  . -  C  ]  =  0.  (4-2) 

o  Nij 

where  ,  C  ,  C^T  denote  the  sub-matrices  of  the  top  row 

of  the  transition  matrix.  The  former  scheme  generally  leads 
to  eigensolution  methods,  whereas  the  latter  generally  leads 
to  root  searching  methods.  An  exception  is  an  application 
of  Laguerre's  root  searching  method  to  a  matrix  in  Kessen- 
burg  form  [20]. 

Laguerre’s  method  and  various  other  root  searching  meth¬ 
ods  are  discussed  in  Section  4.2.  Eigensolution  methods  for 
unsymmetric  matrices  are  presented  in  Section  4.3.  A  survey 
of  eigensolution  methods  for  symmetric  matrices  may  be  found 


[21]  . 


4.2  ROOT  SEARCHING  METHODS 

We  are  required  to  find  X^  ,  i=l,2, . N_(Nt  +  1),  such 

that 

dec  [B±]  =  0  (4-3 

where 


B. 

i 


nt+i  n_ 

[X.  I  -  A  .  C 
1  10 


(4-4) 


is  an  N  xN  polynomial  matrix. 

We  consider  three  techniques  for  obtaining  the  roots  of 
equation  (4-3).  In  the  first  approach,  (4-3)  is  solved  di¬ 
rectly.  This  requires  root  searching  methods  which  utilise 
function  iteration  since  the  explicit  coefficients  of  the 
characteristic  equation  are  not  known.  Muller's  method  [22] 
represents  a  logical  method  for  solution  and  is  discussed  in 
Section  4.2.1.  An  alternate  method  for  obtaining  these 
roots  is  to  use  the  complex  contour  integration  method  of 
Singaraju,  Giri,  and  3aum  [15].  This  technique  is  presented 
in  Section  4.2.2.  The  third  approach  is  to  exploit  polyno¬ 
mial  matrix  reduction  methods  [23]  whereby  the  polynomial 
matrix  (4-4)  is  iteratively  reduced  into  a  triangular  polyn¬ 
omial  matrix.  The  explicit  characteristic  equation  is  then 


e  product  of  the  diagonal  polynomials. 


A  wide  selection 
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of  efficient  polynomial  zero  searching  methods  may  then  be 
used  to  find  the  roots.  A  polynomial  matrix  reduction  meth¬ 
od  is  discussed  in  Section  4.2.3. 

An  application  of  Laguerre's  method  is  presented  in  Sec¬ 
tion  4.2.4.  Since  this  technique  does  not  utilize  either 
form  (4-3)  or  (4-4)  we  consider  it  to  be  independent  of  the 
methods  previously  mentioned.  However,  since  the  method  is 
a  zero  searching  method  it  logically  belongs  within  Section 
4.2. 

4.2.1  Muller ' s  Method 

The  following  is  a  brief  summary  of  the  work  due  to  Mull¬ 
er  [22]  . 

We  are  interested  in  determining  the  values  of  X  which 
satisfy  f(X)=0,  for  some  function  f.  One  begins  the  process 
with  the  values  Xi ,  hi ,  k  ±,  f(xi_1)/  and 

where  X^ ,  hi,  and  ki  are  some  judicious  initial  guesses,  and 
i  is  an  iterative  index;  k^^  is  then  determined  by  the  for- 
mu  1  a 


-2£(.\,)  5 

-L. 


-%i~l 


1  , 


(1-5) 


Si±<V4f  ( V  '5iki  t f  ki"f  <  W  di+f  ( i}  ]  *' : 


where 
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The  process  iteratively  continues  until  some  specified 
criterion  for  acceptance  of  the  root  estimate  is  satisfied 
or  an  upper  limit  on  the  number  of  allowed  iterations  has 
been  reached. 

It  is  interesting  to  note  that  the  convergence  properties 
of  Muller' s  method  have  never  been  proven  for  polynomials 
with  orders  greater  than  two.  Nevertheless  it  is  commonly 
used  on  relatively  large  polynomials  with  excellent  results. 

For  our  purpose,  equation  (4-3)  denotes  the  function  f 
discussed  above.  In  general,  only  one  determinant  evalua¬ 
tion  is  required  for  each  estimate  at  X i-  Excellent  results 
(eight  to  ten  digit  agreement  with  known  solutions)  were  ob¬ 
tained  using  this  technique  for  systems  which  possessed  ap¬ 
proximately  110  roots.  For  higher  order  systems,  however, 
fewer  and  fewer  of  the  predicted  roots  had  any  relation  to 
the  actual  roots.  In  particular,  for  a  system  which  was 
known  to  have  225  roots,  only  4  of  the  predicted  roots  had 
any  relation  to  the  actual  roots.  This  breakdown  is  attri¬ 
buted  to  decreased  separation  in  the  roots  of  large  systems 
(since  by  stability,  all  the  eigenvalues  must  fall  within 
the  unit  circle  in  the  complex  plane),  coupled  with  the  num¬ 
erical  roundoff  errors  associated  with  evaluating  equation 
(4-3) . 
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4.2.2  Contour  Integration 

The  following  brief  summary  is  based  on  the  work  of  S An¬ 
gara  ju,  Giri,  and  3aun  [15]. 

If  a  function  f(X)  is  meromorphic2  in  a  simply-connected 
domain  D  containing  a  Jordan  contour  C,  and  g(X)  is  some 
analytic  function  within  D,  we  may  write  using  the  residue 
theorem 

1  f  f '  (  a  ) 


N 


2^i  Jc  f (X) 


;  ( A ) dA  =  l  g  ( A  ) 

i=l  i 


(4-6) 


where  X0i is  the  i-th  zero  of  f(X)  in  C,  Nq  is  the  total 
number  of  zero's  within  C,  and  the  prime  denotes  differenti¬ 
ation. 

Since  g( X )  is  an  arbitrary  analytic  function,  we  let 
g(X)=X  ,  k=0 , 1 N0 .  The  zero's  of  f(X)  in  C  may  then  be 
obtained  from  the  non-linear  system 
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2  A  meromorphic  function  is  a  function  which  may  be  repre¬ 
sented  as  the  quotient  of  two  entire  functions  and  which 
possesses  poles  only  in  the  finite  complex  plane. 
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where 


\K  y7~-  dA  k  =  0,1,2 . N  (4-8) 

This  rather  unique  technique  has  been  shewn  to  be  quite 
effective  for  obtaining  the  natural  frequencies  from  the 
space-f requency  formulation  of  the  linear,  thin-wire  prob¬ 
lem.  In  this  formulation,  the  function  f(X)  mentioned  above 
corresponds  to  det[Z(sa)],  where  sa  denotes  the  natural  fre¬ 
quencies,  and  Z(«)  denotes  the  moment  method  impedance  ma¬ 
trix.  If,  for  example,  twenty  unknowns  are  used  for  the 
discretization,  then  the  evaluation  of  this  determinate  is 
basically  equivalent  to  the  evaluation  of  a  twentieth  degree 
polynomial.  Numerically,  this  evaluation  should  not  present 
many  complications;  whence,  the  evaluation  of  the  contour 
integral  (4-8)  and  the  subsequent  solution  of  the  non-linear 
system  (4-7)  should  be  numerically  quite  stable. 

The  situation  is  a  bit  more  complicated  for  the  TD-5EM 
formulation  of  this  problem.  The  determinate  of  equation 
(4-4)  now  denotes  the  function  f(X)  above.  For  a  similar 
twenty  unknown  discretization,  f(\)  now  corresponds  to  the 
evaluation  of  approximately  a  four  hundredth  degree  polyno¬ 
mial.  This  represents  a  serious  accuracy  problem  numerical¬ 
ly.  It  was  feared  complications  similar  to  those  observed 


with  Muller's  method  would  develop  with  this  method  when  the 
determinate  of  equation  (4-4)  was  evaluated  prior  to  the 
integration,  and  therefore  the  method  was  not  further  pur¬ 
sued  in  this  context.  By  interchanging  the  determinate  and 
integral  operations,  however,  it  may  he  possible  to  avoid 
the  numerical  errors  associated  with  the  determinate  evalua¬ 
tion.  The  interchange  may  effectively  result  in  a  'numeri¬ 
cal  smoothing'  which  will  give  meaning  to  the  evaluation  of 
the  determinate  even  for  large  scale  problems.  Confirmation 
of  this  conjecture  is  delayed  to  a  future  study. 


4.2.3  Polynomial  Matrix  Reduction 

A  matrix  of  polynomials  may  be  triar.gularized  by  using 
similar  elimination  methods  to  those  associated  with  the  re¬ 
duction  of  standard  matrices  [23].  A  simple  example  is  the 
most  efficient  way  to  describe  the  method. 

Example  4 . 1 

We  consider  the  matrix 


A  I-AC  - C .  = 
o  1 


-  AC  -  C 

"  XC0 

-  C.  i 

°11  X11 

°12 

x12  i 

1 

(4-9) 

-  'C  -  c. 

-  xcn 

-  ci  ! 

°21  *31 

J  />  2 

wnere 
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C  = 
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H 

O 

U 

1 _ 

V 

A, 

V 

!c°n 

A 

=  C1  = 

w  - 

L 

(4-9a) 


By  performing  elementary  operations,  this  matrix  may  be  re¬ 
duced  to  triangular  form  to  yield  the  explicit  characteris¬ 
tic  ecruation 


x"  -  >3  <C0  +co  >  -  x*  tci  +C1  tCo  co  "co  co  > 

22  11  22  111  U12  U21  Uil  22 


+  X  (CQ  C.,  + 

11  “22 


:o  _C0  C1  “C1  C0 
u22  “11  12  21  12  U21 


)  +  (C  c  -c  c  ) 

x22  11  12  21 


(4-10) 


Any  of  a  wide  variety  of  polynomial  zero  searching  methods 
may  now  be  used  to  determine  the  roots. 

This  example  establishes  the  basic  technique.  In  theory 
it  may  be  applied  to  a  matrix  of  arbitrary  size.  Unfortu¬ 
nately,  in  practice  the  method  numerically  breaks  down  due 
to  piling  of  the  coefficients  of  the  eliminated  polynomials 
on  the  diagonal  polynomials.  This  results  in  a  wide  dynamic 
range  in  the  diagonal  coefficients  which  causes  simultaneous 
overspill  and  undersoil!.  This  was  observed  for  systems 
with  only  56  roots.  A  sophisticated  machine  based  scaling 
system  [24]  could  have  partially  controlled  this  dynamic 
range  difficulty;  however,  it  was  feared  it  would  simply 
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postpone  the  breakdown  to  a  slightly  larger  system  and  was 
therefore  not  pursued. 

4.2.4  Laguerre ' s  Method 

An  application  of  Laguerre' s  method  suitable  to  the  ei¬ 
genvalue  problem  has  been  developed  by  Parlett  [20].  The 
following  brief  summary  is  based  on  his  work. 

Let  X  be  an  approximation  to  a  root  of  the  polynomial 
p(\),  where  p(X)  is  of  degree  n.  Laguerre' s  method  requires 
p(X),  p ' ( X ) ,  and  p''(X)  (prime  denotes  a  derivative  with  re¬ 
spect  to  the  argument)  to  obtain  a  better  approximation.  3y 
defining 


?'  (X) 
p(X) 


and  s ,,  ( X ) 


(p 1  (X ) )  ~  -  ?(X)p"  (X) 
(p(X))2 


(4-11) 


Parlett  derives 


a  .  =  a.  -  - - — = —  (.-+-1.:) 

1  ^  1  s1+((n-l)(ns2-sp)1/2 

where  the  square  root  which  maximizes  the  absolute  value  of 
the  denominator  is  chosen  and  n  denotes  the  degree  cf  the 
polynomial  p(X)  (for  the  details  of  this  expression  one 
should  refer  to  Parlett  [ 20 ] ) . 
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Several  convergence  properries  have  been  proven  for  La- 
guerre's  method.  A  listing  nay  be  found  in  Kelly  [25].  . 

Although  the  formulas  given  seem  applicable  only  for  po¬ 
lynomials,  they  may  be  used  on  a  matrix  in  Kessenburg  form 
by  using  Hyman's  method  [20]  to  recursively  yield  the  re¬ 
quired  derivatives  directly  from  this  matrix.  Unfortunate¬ 
ly,  $  is  in  block  upper  Hessenburg  form  and  not  standard  up¬ 
per  Hessenburg  form.  To  use  Hyman's  method  a  transformation 
to  standard  upper  Hessenburg  form  (Section  4.3.1)  would  be 
required.  Such  a  transformation  generally  destroys  the 
sparse  properties  of  a  matrix,  thereby  making  use  of  Laguer- 
re ' s  method  in  this  context  unfeasible  for  determining  the 
required  eigenvalues. 

In  summary,  the  direct  root  searching  methods  tested 
which  exploit  the  form  of  equation  (4-3)  (i.e.,  Muller's 
method,  and  the  polynomial  matrix  method)  were  found  to  be 
useful  only  for  relatively  small  systems  due  to  root  crowd¬ 
ing  and  errors  associated  with  the  numerical  process  (in 
particular,  the  determinant  evaluation) .  The  contour  inte¬ 
gration  technique  may  prove  useful  for  large  systems  if  the 
integration  can  numerically  smooth  the  required  determinate; 
the  feasiblity  of  this  requires  further  consideration.  ?ar- 
lett's  application  of  Laguerre's  method  is  an  excellent  one 
for  solving  large  sparse  matrices  in  Hessenburg  form.  To 


52 


Typically,  the  tri angular! cation  process  requires  several 
seeps  since  only  a  few  elements  of  P  are  operated  on  with 
each  similarity  transformation. 

The  most  common  triangularization  routines  available  are 
the  LR  [27]  and  QR  [23]  algorithms.  Both  of  these  methods 
are  relatively  inefficient  for  fully  populated,  or  dense  ma¬ 
trices  (few  zero  elements).  However,  when  a  matrix  exhibits 
a  certain  pattern  of  zero's  they  become  quite  efficient. 

The  desirable  pattern  of  zero  elements  for  general  matrices 
is  that  which  is  associated  with  an  upper  Hessenburg  matrix. 
To  effectively  use  the  QR  or  LR  algorithms,  then,  one  must 
initially  transform  the  general  matrix  of  interest  to  upper 


Hessenburg  form.  This 


reduction  is  accomplished  by  using 


either  orthogonal  transformations  which  require  approximate¬ 
ly  5n3/3  multiplications,  or  elementary  stabilized  transfor¬ 
mations  which  require  approximately  5n3/5  multiplications 


4. 3. 1.1  The  LR  Transformation 

Let  P,^  be  the  matrix  obtained  from  the  (k-l)-th  transfor¬ 
mation.  ?k  may  be  factored,  or  decomposed,  as 


?,  =  L.  R, 
*  it  V. 


(4-ib) 


I 
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where  L,  is  a  lower  triangular  matrix  with  unit  diagonal 
elements  and  Rk  is  an  upper  triangular  matrix.  The  updated 
iterate  is  then  obtained  by  forming  (k=l,2,...) 

\  he-  (4-15) 

Combining  (4-14)  and  (4-15)  yields 

fk+i  -  ?k  Li  <4-l6> 


which  establishes  the  similarity  transform  form  and  hence 
the  preservation  of  the  eigenvalues  at  each  iteration.  The 
(k+l)-th  transformation  may  be  stepped  back  to  the  original 
matrix  P  by  writing 


Lk-1  V 


(4-17) 


It  can  be  proven  that  the  eigenvalues  of  smallest  moduli 
tend  to  converge  first  and  that  the  rate  of  convergence  is 
dependent  on  the  ratio  of  the  moduli  of  neighboring  eigenva¬ 
lues  [27],  3y  introducing  origin  shifts  the  convergence  rate 
can  be  improved  ( 19 ] . 

Approximately  2  to  3  LR  transformations  are  required  per 
eigenvalue.  As  eigenvalues  are  found,  the  amount  of  re¬ 
quired  computation  steadily  decreases  due  to  smaller  matric¬ 
es  which  must  be  operated  on.  Approximately  n2  multiplica- 
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cions  are  required  for  Che  firsc  few  transformations,  bur 
only  3n3  to  4n3  multiplications  (including  tne  upper  Hessen- 
burg  transformation)  are  required  to  obtain  a  full  set  of  n 
eigenvalues . 

A  thorough  theoretical  discussion  of  the  LR  algorithm  may 
be  found  in  [19];  computational  aspects  may  be  found  in 
[30]  . 

The  use  of  an  orthogonal  factorization  introduces  favora¬ 
ble  stability  and  accuracy  properties  throughout  the  entire 
triangularizaticn  process  [19].  Such  a  decomposition  leads 
to  the  QR  transformation. 

4 . 3 . 1 . 2  The  QR  Trans format! on 

When  the  lower  triangular  matrix  used  in  the  LR  algorithm 
is  replaced  by  an  orthogonal  matrix  we  obtain  che  most  basic 
form  of  the  QR  algorithm.  Let-cing  Q,.  be  the  k-th  iterate 
orthogonal  matrix  we  may  write  the  basic  steps  as 


P,  -  Q,  8. 

k  k  K. 


(4-13) 


and 


W  \  Qk  =  V  ?k\ 


or 


£  X 


(--19) 


I 


The  convergence  properries  of  the  QR  method  are  similar 
to  the  LR  in  that  convergence  is  toward  the  eigenvalues  of 
least  moduli,  and  approximately  1  to  2  double  QR  transforma¬ 
tions  [19]  are  required  per  eigenvalue.  The  orthogonal  de¬ 
composition  does  require  more  computation,  however.  Approx¬ 
imately  5n2  multiplications  are  required  for  the  first  few 
double  transformations;  the  entire  process  requires  approxi¬ 
mately  4n3  multiplications  (including  the  upper  Hessenburg 
transformation) . 

Although  slightly  more  computation  is  required,  the  QR 
method  is  preferred  over  the  LR  due  to  its  superior  stabili¬ 
ty  and  accuracy  properties  for  obtaining  both  single  and 


multitie 


:al  and  complex  eigenvalues. 


. s  tne  case 


when  the  original  matrix  is  real.  A  version  of  both  of 
these  algorithms  exists  for  complex  matrices  [19,3C].  In 
practice,  the  complex  LR  algorithm  has  been  preferred  to  the 
complex  QR  since  it  is  somewhat  simplier  in  content  but  com¬ 
parable  in  stability  and  accuracy. 

The  LR  and  QR  algorithms  represent  the  most  accurate  and 
efficient  methods  available  for  obtaining  a  full  set  of  ei¬ 
genvalues  from  a  dense  matrix  which  has  been  transformed  to 
upper  Hessenburg  form.  A  full  set  of  corresponding  eigen¬ 
vectors  may  be  obtained  by  accumulating  the  transformations 
used  in  the  LR  or  QR  reductions  (this  increases  the  number 
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of  required  multiplications  by  approximately  a  factor  of 
two);  a  partial  set  may  be  obtained  by  using  an  inverse  it¬ 
eration  method  [19].  To  use  the  LR  or  QR  methods,  however, 
matrices  must  generally  be  stored  in  full  storage  mode  since 
similarity  transformations  typically  destroy  any  sparse 
properties.  Hence,  these  methods  are  only  useful  on  matric¬ 
es  with  orders  less  than  a  few  hundred.  The  size  of  the 
transition  matrix  generated  by  TD-SEM  can  easily  be  on  the 
order  of  thousands  and  therefore  excludes  itself  from  full 
eigensolution  by  these  techniques.  Only  partial  eigensolu- 
tion  by  iterative  eigenvector  methods  remain. 

4.3.2  I cerat ive  Eigenvector  Methods 

Iterative  eigenvector,  or  power,  methods  may  be  used  to 
find  either  a  full,  or  more  commonly,  a  partial  set  of  ei¬ 
genvalues  and  eigenvectors.  Since  only  matrix  multiplica¬ 
tions,  in  general,  are  required  by  these  methods  any  sparse 
properties  of  the  original  matrix  may  be  taken  advantage  of. 
Iterative  eigenvector  methods  may  be  divided  into  two  class¬ 
es:  single  vector  and  multiple  vector  methods.  3oth  methods 
begin  with  an  initial  estimate  or  guess  an  an  eigenvector 
which,  hopefully,  will  iteratively  converge  no  an  acnual  ei- 
genvecnor  of  the  system.  The  corresponding  eigenvalue  is 
consequently  found. 


57 


» 


1 


k 


[« 


k 


m 


t 


4.3.2. 1  Single  Vector  Power  Methods 

The  standard  and  inverse  power  methods  represent  the  two 
single  vector  methods  commonly  used  in  practice.  We  will 
initially  consider  the  standard  method. 

Let  P  be  an  unsymmetric  matrix  of  order  n  with  n  indepen¬ 
dent  eigenvectors.  An  arbitrary  vector  u*  0  ’  (where  (o)  de¬ 
notes  the  iteration  number)  may  be  expressed  as  a  linear 
combination  of  these  vectors,  i.e.. 


-(0) 


Clql  +  C2q2  +  +  Cnqn 


(4-20) 


where  q±/  c i=l , 2 ,...., n,  denote  the  eigenvectors  and  ar¬ 
bitrary  coefficients  respectively.  Postmuitiplir.g  ?  by  u‘  5  ) 
yields 


U(1)=  ?U(0)  =  l  C.  ?a.  =  v  X.C.q, 


(4-21) 


i-1 


where  X  denotes  the  eigenvalue  corresponding  to  q±.  If  the 

eigenvalues  can  be  ordered  as  |  X  .  | > | X J •••••> | \  |,  then 

12  n 

u( 1 >  should  represent  an  approximation  to  the  eigenvector  q, 
corresponding  to  the  dominant  eigenvalue  X:  .  This  approxi¬ 
mation  will  iteratively  improve  by  forming 


-(k)  _k  -(0)  .  k„  - 

u  =  P  u'  =  a ,-,q 


lClql  +  X2C2q2 


+  A  *  C  Q 
n  n  'n 


(4-22) 


k  k 

since  for  a  large  enough  k,  |\x|  >>|X2|  .The  ratio  |X1|/|X2| 
determines  how  quickly  the  scheme  will  converge  to  some  spe¬ 
cified  degree  of  accuracy.  Convergence  will  be  quite  rapid 
if  there  exists  good  separation  between  these  eigenvalues. 
Convergence  will  be  very  slow  if  this  ratio  is  close  to  uni¬ 
ty.  A  variety  of  modifications  exist  to  speed  convergence 
under  various  circumstances  [19]. 

_  The  inverse  power  method  requires  more  computation  than 
the  standard  power  method,  but  allows  one  to  approximate  ei¬ 
genvalues  and  eigenvectors  other  than  the  dominant  ones.  We 
begin  the  method  with  some  scalar  a,  and  some  initial  vector 
u<0),  and  consider  the  iterative  system  (k=0, 1, 2, . . . ) 


^(k+1)  =  (P  .  a!)-1  u<k) 


(4-23) 


where  I  is  an  identity  matrix.  Note  that  this  is  the  stan¬ 
dard  power  method  applied  to  the  matrix  (P-al)'1  which  pos- 
esses  the  eigenvalues,  l/(Xi~a),  i=l,2,...,n.  The  method 

converges  to  l/(X.-a),  where  X.  is  the  eicrenvaiue  closest  to 

J  J 

a . 

From  the  convergence  properties  of  the  standard  power 
method,  we  note  the  following  properties  of  the  inverse 
method.  When  a  is  zero,  convergence  is  toward  the  least  do¬ 
minant  eigenvalue  of  P.  When  the  value  of  a  is  clcse  to  an 
eigenvalue,  convergence  will  be  quite  rapid,  but  when  a  is  a 


poor  eigenvalue  estimate,  convergence  may  be  quite  slow. 
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In  practice,  the  inverse  observed  in  equation  (4-23)  is 
rarely  explicitly  determined.  A  decomposition  or  Gaussian 
elimination  is  used  on  the  matrix  (P-al)  which  requires  ap¬ 
proximately  n3/3  multiplications.  Another  n2  multipiica- 

-(k) 

tions  are  required  to  determine  each  new  u  (eigenvector 
iterate).  The  standard  power  method  does  not  necessitate  a 
decomposition  and  is  therefore  more  efficient  than  the  in¬ 
verse  method  for  the  same  number  of  iterations.  The  number 
of  iterations  required  for  the  inverse  method,  however,  will 
be  significantly  smaller  when  a  good  guess  at  an  eigenvalue 
is  known.  Typically,  the  two  methods  are  used  in  conjunc¬ 
tion  with  one  another.  The  standard  method  determines  a 
good  initial  guess  which  is  then  refined  by  the  inverse 
method. 

Once  a  single  eigenvalue  and  eigenvector  is  known  it  may 
be  filtered  out  of  the  original  matrix  by  either  purifica¬ 
tion  or  deflation  [29].  The  power  methods  above  may  then  be 
used  on  this  filtered  matrix  to  find  another  eigenvalue  and 
eigenvector.  These  may  then  be  filtered  out,  and  so  on. 

The  state  transition  matrix  is  not  only  very  sparse,  but 
possesses  a  full  set  of  eigenvalues  whose  moduli  are  less 
than  unity.  These  two  properties  exclude  the  use  of  single 
vector  power  methods  for  the  following  two  reasons:  first, 
the  purification  or  deflation  processes  which  are  required 


to  find  several  eigenvalues  possess  the  unfortunate  property 
that  they  destroy  sparsity  in  general;  second,  the  ratio  of 
the  moduli  of  neighboring  eigenvalues  is  generally  so  close 
to  unity  that  convergence  is  impractical  in  a  reasonable  am¬ 
ount  of  time.  Experimental  verification  of  these  two  com¬ 
plications  has  established  that  single  vector  power  methods 
do  not  constitute  a  feasible  method  of  solution  for  the 
large  scale  state  transition  matrix. 

4. 3. 2. 2  Multiple  Vector  Power  Methods 

Multiple  vector  methods  are  power  methods  which  iterate 
with  several  vectors  simultaneously.  The  name  simultaneous 
;t=rar;cn  (5.1.)  has  beer,  given  to  these  methods. 

Bauer  [31]  introduced  the  first  S.I.  concept,  called 
' Bi- Iteration' .  This  method  solved  the  algebraic  eigenvalue 
problem  ?u=\u  for  an  arbitrary  matrix  P  of  order  n.  The 
idea  of  the  method  was  to  iterate  with  two  sets  of  vectors, 

u1,u2,  .  .  .  .  ,us  and  q^q,, . ,qs  (s£n)  applied  to  ?  and  ?n 

(H  denotes  Hermitian  transpose)  respectively.  By  maintain¬ 
ing  u^_  and  biorthonormal  where  6^  is  the  Kro- 

necker  delta),  it  can  be  shown  [31]  that  under  certain  con¬ 
ditions  the  u,.  converge  to  the  right  eigenvectors,  and  the  q 


ration,  emphasis 
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We  introduce  the  following  geometrical  argument  to 
justify  this  sequence  of  operations.  Let 


9 


Vk  "  P  Uk-1 


(4-25) 


fl 


where  U, contains  an  approximation  to  a  subset  of  the  ei¬ 
genvectors  of  P.  may  be  decomposed  into  the  sum  of  a 

projected  matrix  and  an  orthogonal  matrix. 


M 


v,  =  v. 

k  tc 


rproj  .  +  vortho. 


M  (it  D  )  +  V°rCh0 ' 

ak-i  V  T  vk 


_  „  T  „  ,  ..ortho. 

=  a.  T.  D.  +  v. 

\  K  r(.  iC 


where  Ek,  Z-  ,  and  Dk  are  to  be  determined.  Combining  (4-25) 
and  (4-26)  yields 


?U.  .  =  E.T,D.  -  V°rth°- 
K.-1  &  &  k  s 


if  and  only  if 
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(4-28) 
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Equation  (4-29)  defines  an  eigenvalue  problem 


for  Dj.  ;  the  eigenvalues  of  Dk  are  given  by  the  matrix  0^, 
while  the  eigenvectors  of  Dk  are  given  by  the  matrix  Afc.  D, 
is  found  by  solving  the  linear  system 


U.k  ,  v.  -  (U?  .  U,  .)  D. 

it— 1  k  k-1  k-1  k 


(4-30) 


The  updated  iterate  is  then 


\  -  \  \  ’  \-x  V  +  vkrth0' 


(4-31) 


where  is  a  normalized  W^.  Note  that  in  the  iteration  cy¬ 
cle  an  approximation  which  ignores  the  contribution  of  the 
second  term  on  the  right  of  equation  (4-31)  is  made.  This 
approximation  is  reasonable  since  the  second  term  tends  to 
zero  as  the  number  of  iterations  becomes  large  (due  to  A.a 
approaching  a  span  of  the  projection  space). 

This  sequence  of  steps  is  appropriate  for  a  right  eigen¬ 
vector  or  projection  scheme.  A  similar  sequence  may  be 
developed  for  a  projection  scheme.  In  the  latter  case 

convergence  is  toward  the  reciprocal  eigenvalues. 

A  few  comments  are  in  order  about  the  sequence  of  steps 
in  the  iteration  cycle.  U0  represents  the  initial  guess  to 
a  dominant  block  of  eigenvectors  which  may  be  with  real  num¬ 
bers.  The  entire  process  proceeds  in  real  arithmetic  until 
the  'interaction  matrix',  D.  ,  possesses  complex  eigenvalues; 


the  process  then  becomes  complex.  When  this  transition  is 
made,  becomes  a  positive  definite  Hermitian  matrix  for 
which  a  modified  Cholesky  decomposition  has  been  developed 
to  minimize  the  operations  required  for  the  solution  of  step 
iv. 

If  U,  contains  s  approximations  to  s  eigenvectors  (sSn, 
where  n  again  is  the  order  of  the  original  matrix  P),  then 
the  matrix  0,K  contains  s  approximations  to  s  eigenvalues. 

It  can  be  shown  [34]  that  the  rate  of  convergence  of  the 
above  sequence  is  dependent  on  the  ratio  ] \ 1 1 / ] X  +  , |  and  not 
on  the  ratio  ( X x I / 1 X2 ]  which  governs  the  single  vector  meth¬ 
ods.  Convergence  is  toward  the  dominant  eigenvalues. 

For  certain  geometries  (Section  5.3.2)  the  eigenvalues 
associated  with  the  state  transition  matrix  exibit  a  well 
defined  series  of  annular  rings.  Within  each  of  these  rings 
the  eigenvalues  are  of  comparable  moduli.  3etween  rings 
there  is  a  significant  jump  in  moduli.  Since  the  conver¬ 
gence  rate  of  the  simultaneous  iteration  scheme  above  is  de¬ 
pendent  on  the  ratio  of  neighboring  blocks  of  eigenvalues 
instead  of  the  ratio  of  neighboring  eigenvalues,  and  since 
the  scheme  requires  only  matrix  multiplications  which  pre¬ 
serve  sparsity,  simultaneous  iteration  represents  a  reason¬ 
able  method  for  obtaining  a  partial  set  of  eigenvalues  and 
eigenvectors  for  $.  Theoretically,  rings  may  be  found  one 
by  one  until  a  desired  set  is  collected. 


A  complication  with  this  method  is  the  number  of 


multiplications  required  per  iteration.  For  a  real  arith¬ 
metic  cycle  approximately  [29] 


,5  2  31  3 
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c  0 


(4-32) 
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multiplications  are  required.  Here,  c  denotes  the  average 
number  of  nonzero  elements  in  each  row  of  the  matrix  whose 
eigensolution  is  sought,  n  denotes  the  order  of  this  matrix 
and  s  denotes  the  number  of  eigenvalues  desired.  For  a  com 
plex  cycle  this  number  increases  approximately  by  a  factor 
of  four. 

As  the  order  of  the  matrix  whose  eigenvalues  are  scucht 
increases,  the  separation  between  neighboring  blocks  of  ei¬ 
genvalues  typically  decreases,  and  therefore,  more  itera¬ 
tions  are  required  to  obtain  a  specified  degree  of  accuracy 
It  can  be  shown  that  the  approximate  number  of  required  it¬ 
erations  (k)  to  obtain  a  certain  degree  of  accuracy  is 

k  “  lne/lnj,\1/Xs+1i  (4-33) 

where  s  denotes  rhe  block  size,  and  e  denotes  the  accuracy 
desired, . i . e . ,  0 . CC01 , 0 . 001,  ere. .  Figure  5.9  (Section 
5.3.2)  illustrates  size  of  the  state  transition  matrix  vs. 
execution  time  required  on  an  IBM  3032  computer  for  a  parti 


cular  subset  of  eigenvalues.  Figure  5.10  illustrates  sto¬ 
rage  requirements  for  both  the  QR  and  S.I.  methods. 

A  significant  reduction  in  the  number  of  multiplications 
required  per  iteration  may  be  accomplished  by  taking  advan¬ 
tage  of  the  form  of  the  eigenvectors  associated  with  the 
transition  matrix.  For  simplicity,  equation  (3-32)  is  re¬ 
peated  here  as 


— T  r  “T— T 

K  *  LZ,.  V  ,  Z  V  . .  V  ] 

Cl  a  a  a  a  ’  *  aJ 


(4-34) 


where  z  denotes  a  particular  eigenvalue,  and  v  denotes  a 

'*•  Cl 

vector  spatially  describing  the  natural  mode.  Since  the  en¬ 
tire  state  vector  X  can  be  constructed  once  these  two  Dar- 

a 

ameters  are  known,  it  should  be  possible  to  carry  out  the 
simultaneous  iteration  cycle  with  only  a  subset  of  the  state 
vector.  Under  certain  conditions  this  is  possible  and  has 
been  given  the  name  'sub'  iteration. 

The  sub  iteration  modification  may  only  be  used  after  the 
proper  form  of  the  entire  state  vector  has  appeared.  This 
can  typically  take  several  iterations  since  the  vectors  are 
constructed  from  the  top  down.  A  sophisticated  computer 
program  which  incorporates  this  idea  has  shown  that  the  de¬ 
sired  form  cannot  be  forced  in  general,  it  must  appear  na¬ 
turally.  Once  the  form  appears,  however,  the  use  of  sub  it¬ 
eration  for  the  remaining  iterations  will  drastically  reduce 
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the  execution  time.  A  listing  of  this  routine  (which  was 
designed  for  the  solution  of  the  thin-wire  problem  (Chapter 
5))  may  be  found  in  Appendix  A.  Random  sparse  storage  tech¬ 
niques  [29]  were  used  in  this  algorithm  to  store  the  state 
transition  matrix. 

Of  all  methods  tested,  simultaneous  iteration  with  the 
sub  iteration  modification  represents  the  only  feasible 
method  for  the  partial  eigensolution  of  the  large  scale 
state  transition  matrix.  The  amount  of  computation  time  re¬ 
quired  using  this  scheme  can  admittedly  become  excessive  as 
the  system  size  becomes  very  great.  However,  it  should  be 
realized  that  even  if  it  were  possible  to  place  these  large 
matrices  in  the  high  speed  store  of  the  largest  computers  so 
that  the  efficient  QE  or  LR  algorithms  could  be  used  to  find 
both  the  eigenvalues  and  eigenvectors,  the  amount  of  compu¬ 
tation  time  would  be  at  least  as  great  if  not  greater  than 
the  time  required  by  S.I..  This  may  be  established  simply 
by  considering  the  total  number  of  operations  each  scheme 
requires.  The  time  required  by  the  QR  and  LR  methods  would 
be  for  a  full  set  of  eigenparameters  as  opposed  to  the  par¬ 
tial  set  given  by  S.I.,  however. 


Chapter  V 

WIRE- STRUCTURE  ANALYSIS 


5.1  INTRODUCTION 

Wire  structures  may  be  analyzed  either  as  scatterers  or 
antennas.  When  an  incident  wave  propagating  in  space  ex¬ 
cites  a  response  on  the  wire,  we  consider  the  structure  to 
be  a  scatterer.  When  the  wire  is  excited  from  small  regions 
on  the  structure  itself,  it  is  considered  an  antenna. 

In  this  chapter,  we  will  consider  borh  frequency-  and 
time-  domain  numerical  solution  techniques.  The  discussions 
will  be  restricted  to  thin,  perfect  conducting  wires  situat¬ 
ed  in  one-dimensional  or  linear  geometries.  Frequency-domain 
methods  are  presented  in  Section  5.2.  By  applying  the 
widely  used  moment  method  [14]  to  the  EFIE  for  thin,  perfect 
conducting  surfaces,  the  current  and  electric  field  distri¬ 
butions  along  the  antenna  are  obtained.  A  discussion  of  ex¬ 
citation  models  is  presented  in  Section  5.2.1;  and  a  discus¬ 
sion  of  the  applicability  of  the  thin-wire  kernel 
approximation  is  presented  in  Section  5.2.2.  In  Section 
5.3,  time-domain  techniques  are  discussed.  3v  using  finite 
difference  approximations,  we  obtain  the  transient  current 
distribution  associated  with  the  antenna  mode  of  operation; 
and  by  applying  the  TD-SEM  method,  we  present  the  SEM  pole 
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distribution  associated  with  the  scattering  mode.  The  ef¬ 
fect  of  varying  the  time-sampling  distance  used  in  the  fi¬ 
nite  difference  approximation  on  the  pole  locations  is  pre¬ 
sented  in  Section  5. 3. 2.1.  In  Section  5.3.3  ,  the  effect  of 
segment  to  segment  (unknown  to  unknown)  coupling  through  the 
kernel  of  the  EFIE  on  the  pole  locations  is  discussed. 

\ 

Extensive  studies  on  approximate  analytic  solution  tech¬ 
niques  of  linear  wire  structures  have  been  given  by  King 
[35].  Numerical  solution  techniques  have  been  given  by  Har¬ 
rington  [14 1,  Poggio  and  Mayes  [36],  Thiele  [37],  and  Mittra 
[38].  These  provide  excellent  discussions  of  methods  for 
obtaining  most  any  desired  antenna  or  scattering  parameter 
ocher  chan  the  pole  structure.  A  complete  discussion  and 
several  references  for  obtaining  the  natural  responses  of 
objects  in  both  the  frequency-  and  time-domains  may  be  found 
in  Chanter  3 . 


5.2  SPACE -FREQUENCY  TECHNIQUES 

The  appropriate  expression  which  describes  thin,  perfect 
conducting  structures  is  (2-lla).  For  one-dimensional  prob¬ 
lems  with  geometries  similar  to  the  geometry  depicted  in 


Figure  5.1,  this  expression  may  be  written  in  terms  of  an 
unknown  current,  I,,  and  magnetic  vector  potential.  A,,  as 


-E““'"(a,z  ,y)-- 


!  k"+— ~  ;  A2(a,z;-0 


(5-1) 


We  let  the  structure  of  interest  be  uniformly  divided 
into  (N+l)  pieces  each  of  length  A.  The  current  distribu¬ 
tion  may  then  be  approximated  as 


i’(z’)  -Vis 
L.  n  n 
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denotes  the  piecewise  sinusoids,  and  the  constants  I  are  to 
be  determined.  3y  substituting  this  approximation  into 
equation  (5-2),  we  may  write  a  discrete  version  of  expres¬ 
sion  (5-1),  after  several  straightforward  manipulations,  as 


-Elnc  =  Til  (s  ) 
z  n  od  n 

n*l 


(5-6) 


where  the  operator  Lop(sn)  is  given  by 
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Here,  the  arguments  of  the  function  G  replace  (z-z')  in 
either  equation  (5- lb)  or  (5-4a)  depending  on  whether  the 


exact  or  aoorcximate  kernel  is  used. 
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The  matrix  with  components  <w  , L  (S  )>  is  known  as  an  imoe- 

^  m  op  '  a7 

~  i  r\  r 

dance  matrix  2,  and  the  vector  with  components  <w  , -a  >  is 

a  2 

known  as  a  voltage  vector  V.  Using  this  notation,  the  cur¬ 


rent  coef f icier.ts  mav  be  found  bv 


I  =  z_1V 
n 


(5-10) 


assuming  a  non-singular  impedance  matrix. 

All  frequency-domain  results  in  this  chapter  will  be 
based  on  this  model. 


.2.1  axcitarion  Models 

The  impedance  matrix  is  invariant  to  me  form  of  excita- 


moae,  the  excitation  may  be  modeled  in  a  varie*  t  o-  •  ays 
The  simplest  is  the  delta  function  [37]  gener-  ..„i 


oltace  vector  reflects  ohm.  In  the  antenna 


which  introduces  a  single  entry  of  unity  into  the  voltage 
vector;  its  position  in  the  vector  corresponds  to  the  posi¬ 
tion  of  the  excitation  on  the  antenna.  An  alternate,  and 
more  rigorous,  method  of  excitation  modeling  is  the  use  of 
equivalent  magnetic  sources.  This  leads  to  magnetic  frill 
[37]  generators  and  belt  [40]  generators.  The  voltage  vec¬ 
tor  reflects  this  form  of  excitation  in  the  form  of  a  dis¬ 
tribution  with  unit  area.  Results  induced  by  these  diffe¬ 
rent  models  are  nearly  identical  everywhere  over  the 
structure  except  directly  at  the  point  where  the  excitation 
is  originating.  This  point  is  somewhat  critical  to  the  over 
all  analysis;  however,  since  impedence  and  admittance  char¬ 
acteristics  are  dependent  on  the  magnitude  of  the  current  at 
this  point.  It  has  been  found  [40]  that  the  magnetic  source 
generator  models  yield  impedance  and  admittance  values  in 
better  agreement  with  experimental  measurements  than  delta 
generator  models.  Figure  5.2  compares  a  delta  generator 
with  a  magnetic  belt  generator  for  a  half-wavelength  antenna 
using  the  exact  kernel.  Numerically,  the  results  differ 
only  at  and  near  the  feed  point;  graphically  the  difference 


is  indistinguishable . 
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5.2.2  Thin- Wire  Kernel  Approximation 

Equation  (5-4)  denotes  the  thin-wire  kernel  approximation 
which  is  analytically  valid  for  wires  that  are  only  a  small 
fraction  of  a  wavelength  in  diameter.  A  numerical  complica¬ 
tion  can  arise  using  this  kernel  even  when  the  wire  is  elec¬ 
trically  and  hence  the  approximation  analytically  valid. 

When  the  sampling  distance  between  unknowns  is  on  the  order 
of  the  radius  of  the  structure,  non-physical  current  oscil¬ 
lations  may  appear  both  at  the  end  points  and  the  feed 
point.  Use  of  the  exact  kernel  in  such  cases  avoids  this 
complication.  Figure  5.3  compares  the  effect  of  the  exact 
and  approximate  kernels  on  the  current  distribution  of  a  1 
meter  antenna.  A  belt  generator  was  used  for  the  excita¬ 
tion. 


5.2.3  Electric  Field  Distribution 

Once  the  current  coefficients  are  obtained  from  equation 
(5-10),  the  scattered  tangential  electric  field  distribution 
may  be  found  by  evaluating 
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Note  that  the  approximate  thin-wire  kernel  has  been  used. 
Figure  5.4  shows  E®  ( 2 )  for  a  1  meter  antenna.  A  delta  gen¬ 
erator  at  the  origin  was  used  for  the  excitation. 


5.3  SPACE -TIME  TECHNIQUES 

For  one-dimensional  structures  corresponding  to  Figure 
5.1,  expression  (2-12a)  may  be  written  as 
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where  the  approximate  thin-wire  kernel  (equation  (5-4))  has 
been  used,  and  the  current  density  appearing  in  expression 
(2- 12a)  has  been  replaced  by  the  total  current. 

The  discretization  of  this  expression  will  follow  the 
technique  developed  in  Section  3.3.1.  We  assume  (from  ex¬ 
pansion  (3-17))  that  the  current  distribution  may  be  approx¬ 
imated  as 
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where  N  denotes  the  number  of  unknowns  used  in  the  discreti¬ 
zation.  With  this  expansion,  the  discrete  vector  potential 
is  given  by 
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Note  that  the  time  pulse  has  been  pulled  out  of  the  integ¬ 
ral.  This  is  valid  only  for  very  thin,  linear  structures. 

By  letting  Ga  denote  the  integral  appearing  in  (5-14), 
and  picking  p'=p-|n-n'|  (note  that  this  assumes  the  choice 
cAt=A),  we  may  write  the  vector  potential  as 
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where 
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By  transforming  the  continuous  differential  operators  ap¬ 
pearing  in  expression  (5-12)  to  central  finite  difference 
operators,  we  obtain  the  following  explicit  scheme  for  the 


unknown  current  coefficients: 
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where  G„  is  the  self  patch  kernel,  and  I_1  ,  IN+1  ,  IN+2  ,  Ig 
are  defined  to  be  zero. 

We  now  consider  the  choice  cAt=A  by  performing  a  Fourier 
stability  analysis  (Section  3.3.1)  on  the  homogeneous  dif¬ 
ferential  equation  for  the  vector  potential.  The  vector  po 
tential  is  analyzed  for  simplicity;  the  stability  require¬ 
ment  for  the  current  equation  is  anticipated  to  be  similar 
due  to  integral  relation  between  the  two.  The  following 
discussion  outlines  the  technique. 

We  begin  with  the  one-dimensionai  wave  equation  for  the 
vector  potential 
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which  in  explicit  differenced  form  becomes 
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where  r2=(cAt/A)2,  and  Az  =A  z(  nA ; pAt )  .  Discrete  separa¬ 

te,? 

tion  of  variable  techniques  then  leads  us  to  a  solution  of 
the  form 

A  =  exp  (jctnA)  (5-19) 

Zn,p 

where  w=exp { sAt }  (s  is  arbitrary  and  may  be  complex),  and  a 
is  an  arbitrary  real  number.  Clearly,  for  this  solution  to 
remain  bounded  for  all  p,  the  magnitude  of  u  must  be  bounded 
above  by  unity.  By  substituting  this  general  solution  into 
the  basic  equation  (5-18),  we  obtain  the  following  quadratic 
equation  in  u>: 

oj2  -  2  (l-2r2sin2(aA/2))u  +  1  =  0  (5-20) 

From  this  equation,  the  magnitude  of  u>  is  less  or  equal  to 
unity  for  r£l.  Therefore,  the  choice  cAt=A  may  yield  a  sta¬ 
ble  solution  for  the  current  coefficient  difference  equation 
(equation  (5-16)).  Certainty  cannot  be  ascertained  for  the 
following  two  reasons:  first,  any  appropriate  boundary  con¬ 
ditions  have  not  been  included  in  the  analysis;  and  second, 
a  difference  equation  for  the  vector  potential  has  been  con¬ 
sidered  rather  than  a  difference  equation  for  the  current 
coefficients.  The  matrix  stability  method  (Section  3.3.1) 
must  be  applied  to  the  current  difference  scheme  for  this 
criterion  to  be  numerically  rigorous. 
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5.3.1  Transient  Current  Response 

Figure  5.5  depicts  the  transient  current  response  ob¬ 
tained  from  equation  (5-16)  on  one  of  two  feed  segments  for 
a  1  meter  dipole  antenna.  The  mathematical  representation 
of  the  excitation  used  on  each  of  these  segements  is  given 


T  exp(-s2^-cn 


(5-21) 


where  Ea  is  the  free-soace  imoedance  120tt,  t  is  the  time 

*  max 

when  the  magnitude  of  the  pulse  peaks,  and  g  is  the  compres¬ 
sion  factor  of  the  pulse.  The  parameter  t^aY  was  chosen  to 
be  0.5  light  meters  (LM),  and  g  was  chosen  so  that  the  mag¬ 
nitude  of  the  pulse  at  t=C  and  t=l  L24  was  C .  0001/(  60tt) 
volts/meter . 


5.3.2  TD-SZM  Pole  Pi szribution 

The  TD-SEM  technique  developed  in  Section  3.3.1  may  be 
applied  to  the  difference  scheme  of  equation  (5-16).  The 
eigenvalues  of  the  state  transition  matrix  i  may  be  found  by 
the  methods  given  in  Chapter  4.  The  order  of  I  for  this 
problem  is  N(N+1),  where  N  represents  the  number  of  unk¬ 
nowns.  For  NS13  full  eigensolution  by  the  QR  transformation 
(Section  4.3. 1.2)  is  recommended.  For  larger  values  of  N, 
partial  eigensolution  by  simultaneous  iteration  (Section 
4. 3. 2. 2)  is  recommended.  Figure  5.6  shows  the  pole  distri- 
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gments  of  a  1  meter  antenna  with  a  radius  of 
00754  meters.  The  number  of  unknowns  was  22. 
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bution  obtained  by  TD-SEM  (QR  solution)  for  a  1  meter  scat- 
terer  with  a  radius  of  0.005  meters  discretized  with  18  unk¬ 
nowns  such  that  A=1/(N+1 )=1/19 .  The  results  have  been  sca¬ 
led  by  1/(ttc),  and  a  comparison  with  the  frequency-domain 
results  of  Singaraju,  Giri  and  Baum  [15]  has  been  made.  The 
results  shown  reflect  only  the  second  quadrant  of  poles 
since  the  complex  conjugates  may  be  obtained  by  symmetry. 

A  few  remarks  on  the  structure  of  the  eigenvalues  associ¬ 
ated  with  this  problem  for  the  choice  cAt=A  should  be  made. 
From  Suability,  all  eigenvalues  must  fall  within  the  unit 
circle  in  the  complex  plane.  Complex  eigenvalues  appear  in 
groups  of  four  which  graphically  define  a  square.  Several 
groups  of  nearly  identical  moduli  form  annular  rings.  Real 
or  purely  imagine' eigenvalues  appear  in  pairs  of  equal  mo¬ 
duli  but  differing  sign.  The  eigenvalues  of  interest  are 
those  with  positive  real  components.  Those  with  negative 
real  components  are  conjectured  to  correspond  to  false  poles 
which  have  no  true  physical  meaning  (this  conjecture  is  dis¬ 
cussed  in  Section  5. 3. 2.1).  The  layering  structure  of  the 
poles  corresponds  to  the  annular  ring  structure  of  the  ei¬ 
genvalues.  For  an  even  number  of  unknowns,  2N  of  the  eigen¬ 
values  in  the  outermost  annular  ring  correspond  to  the  sig¬ 
nificant  poles  contained  within  the  first  layer,  while  2N-2 
eigenvalues  correspond  to  the  poles  of  the  second  layer. 


I 


|  OH.)/ j  S}U1I 


39 


Further  layers  which  are  well  defined  have  successively  4 
fewer  eigenvalues  than  the  preceeding  layer.  For  an  odd 
number  of  unknowns  a  factor  of  2  must  be  added  to  the  num¬ 
bers  above.  Figure  5.7  depicts  the  eigenvalue  structure  for 
18  unknowns.  The  quadrant  which  is  blocked  off  denotes  the 
eigenvalues  of  interest. 

Figure  5.8  shows  the  first  layer  pole  distribution  for  32 
unknowns  as  found  by  simultaneous  iteration  techniques  with 
the  sub  iteration  modification.  It  should  be  noted  that.  a. 
simultaneous  iteration  algorithm  may  be 

used  to  find  more  layers  than  just  the  first.  Additional 
layers  may  be  found  on  a  single  computer  execution  or  multi¬ 
ple,  independent  executions. 

Figure  5.9  compares  execution  time  requirements  (on  an 
I3M  3032  computer  using  FORTRAN  H  EXTENDED  (OFT=2))  with  the 
order  of  the  state  transition  matrix  for:  full  eigensolution 
by  the  QR  algorithm  (eigenvalues  only),  QR  algorithm  (both 
eigenvalues  and  eigenvectors),  and  first  layer  eigenvalues 
and  eigenvectors  by  simultaneous  iteration  (assuming  complex 
iteration  cycles)  with  and  without  the  sub  iteration  modifi¬ 
cation.  IBM  double  precision  was  found  to  be  necessary  when 
the  QR  algorithm  was  used;  IBM  single  precision  was  suffi¬ 
cient  for  the  simultaneous  iteration  method.  The  curves  in 
Figure  5.9  reflect  these  precision  requirements. 


Figure  5.10  compares  storage  requirements  with  the  tran¬ 
sition  matrix  order  for  the  cases  above.  The  IBM  douple 
precision  storage  requirement  for  the  QR  algorithm,  eigenva¬ 
lues  only,  is  8n2  (where  t\  denotes  the  order  of  4).  This 
requirement  becomes  8ii2  +  16ti2  when  both  the  eigenvalues  and 
all  eigenvectors  are  desired.  The  single  precision  storage 
requirement  of  simultaneous  iteration  is  4(4k2+(3+2N)k 
+2nk+3( (2N-1)N) ) ,  where  k  denotes  the  number  of  eigenvectors 
sought,  and  N  denotes  the  number  of  unknowns  (for  Figure 
5.10,  k=2N  was  chosen;  this  defines  the  first  layer). 

5.3.2. 1  Effect  of  Varying  the  Time  Sampling  Distance 
The  pole  structure  presented  in  Figures  5.6  and  5.7  is 
for  the  choice  cAt=A.  When  cAt  is  chosen  to  be  less  than  A, 
complications  arise.  Since  each  time  step  is  smaller,  more 
time  steps  are  required  for  a  wave  to  travel  the  length  of 
the  structure.  This  causes  the  size  of  the  state  transition 
matrix  to  increase,  and  thereby  possess  a  larger  number  of 
eigenvalues.  The  relation  of  these  extra  eigenvalues  to  the 
true  poles  of  the  system  is  of  interest.  It  has  been  found 
that  the  additional  eigenvalues  add  additional  'poles'  which 
may  be  divided  into  two  types.  The  first  type  are  complex 
poles  with  imaginary  components  which  are  conjectured  to  be 
of  greater  magnitude  than  the  'true  system  poles';  the  sec- 


ond  type  are  purely  real  poles  which  are  conjectured  zo  ex¬ 
tend  further  down  the  real  axis  rhan  the  true  poles.  Both 
of  these  conjectures  and  an  elimination  technique  for  the 
'spurious'  poles  are  discussed  in  the  following  paragraphs. 

The  angle  associated  with  the  polar  representation  of  ei¬ 
genvalues  which  possess  positive  real  and  imaginary  compo¬ 
nents  is  at  most  it/2  radians  (by  symmetry,  we  may  similarly 
consider  the  conjugate  of  these  eigenvalues).  The  imaginary 
component  of  the  poles  corresponding  to  these  eigenvalues, 
then,  have  a  magnitude  which  is  at  most  l/(2cAt)  radians 
(recalling  that  the  poles  are  scaled  by  1/(tc)).  When  the 
stability  condition  cAt=A  is  valid,  this  maximum  becomes 
1/ ( 2A )  radians.  Since  all  thin- wire  pole  distributions  ob¬ 
tained  by  frequency-domain  techniques  are  bounded  by  1/(2A) 
radians,  it  is  sensible  to  bound  time-domain  methods  by  that 
value  as  well.  Hence,  we  conjecture  that  the  true  system 
poles  obtained  by  TD-SEM  are  bounded  on  the  imaginary  axis 
by  1/(2A)  for  arbitrary  choice  of  cAt  (note  that  this  res¬ 
triction  eliminates  the  poles  corresponding  to  eigenvalues 
with  negative  real  components).  This  technique  generates  a 
simple  criterion  from  which  spurious  poles  off  the  real  axis 
may  be  eliminated.  We  consider,  next,  the  spurious  poles 
which  are  situated  on  the  real  axis. 


The  second  type  of  additional  poles  are  not  as  simple  to 
emove.  No  well  defined  method  has  been  developed.  These 


poles  do  not  disturb  the  well  defined  layering  structure 
since  they  appear  past  it.  TD-SEM  generates  true  poles  past 
the  layering  structure  (Figure  5.6)  however,  and  therefore 
an  ambiguity  exists.  One  possible  restriction  would  be  to 
retain  only  those  poles  within  the  range  0  to  approximately 
-1/(2A)  on  the  real  axis.  Although  this  criterion  is  admit¬ 
tedly  strict,  it  may  be  effectively  used  for  arbitrary 
choices  of  cAt. 

To  conclude  this  section,  it  should  be  noted  that  the 
left  half  plane  pole  structure  obtained  from  an  unstable  al¬ 
gorithm  has  very  little  similarity  with  the  pole  structure 
obtained  from  a  stable  algorithm.  Hence,  one  should  be  cer¬ 
tain  that  a  particular  scheme  is  stable  if  the  results  ob¬ 
tained  are  to  be  considered  meaningful  and  accurate. 

5.3.3  Pole  Shift  by  Kernel  Decoupling 

It  was  suspected  a  priori  that  a  relation  may  exist  between 

the  sub-  matrices  of  the  transition  matrix  and  the  par¬ 
ticular  layering  structure  of  the  poles.  In  other  words, 
the  first  few  sub-matrices  may  contribute  the  poles  of  the 
first  layer,  the  next  few  the  poles  of  the  second  layer,  and 


30  on. 


To  test  this  conjucture,  the  effect  of  zeroing  the 


kernel  G.  , .  of  equation  (5-16)  for  various  values  of  the 
I  n-n  I 

difference  |n-n' |  was  considered.  This  effectively  removes 
element  to  element  coupling  and  thereby  zeros  sub-matrices. 
Figure  5.11  depicts  the  movement  of  the  pole  distribution 
for  a  1  meter  antenna  with  a  radius  of  0.005  meters  discre¬ 
tized  with  10  unknowns  (N=10).  All  layers  tend  to  shift  to¬ 
ward  the  imaginary  axis  (approaching  a  transmission  line)  as 
elements  are  zeroed  from  the  extreme  outer  sub-matrices  in¬ 
ward.  Unfortunately,  no  relation  between  the  block  sub-ma¬ 
trices  and  particular  layers  can  be  inferred  from  this  form 
of  movement.  Specifically,  the  results  shown  in  Figure  5.11 
may  be  interpreted  as  follows:  full  set  denotes  no  decou¬ 
pling,  the  level  1  set  has  the  kernel  evaluated  at 
|n-n' |=N+1  set  equal  to  zero,  the  level  2  set  has  the  evalu¬ 
ation  at  | n-n' |=N+1  and  N  equal  to  zero,  the  level  3  set  has 
rhe  evaluation  at  | n-n' |=N+1,  N,  and  N-l  equal  to  zero,  and 
the  level  7  set  is  similarly  defined. 


Chapter  VI 


TRANSIENT  ANALYSIS  OF  THIN,  PERFECT  CONDUCTING 
RECTANGULAR  PLATES 

6.1  INTRODUCTION 

The  rectangular  plate  falls  within  the  class  of  struc¬ 
tures  known  as  open  structures  with  edges.  Structure.-  ,.'_th 
edges,  and  in  particular  corners,  are  difficult  to  am  ze 
due  to  the  singular  behavior  of  the  current  component  mal¬ 
lei  to  an  edge,  and  the  ambiguity  of  the  current  magni  -ue 
in  a  corner.  These  complications  have  generally  restricted 
exact  analytic  solutions  to  infinite  half-plane  problems. 

In  particular,  we  cite  an  exact,  frequency-domain  solution 
for  the  current  density  generated  on  an  infinite  half-plane 
by  an  edge  on  incident  plane  wave  which  may  be  found  in  Born 
and  Wolf  [41];  the  transformed,  time-domain  result,  may  sub¬ 
sequently  be  found  in  Davis,  et  al  [42]. 

Since  the  realm  of  problems  which  are  solvable  by  analyt¬ 
ic  methods  is  quite  narrow,  interest  has  turned  -toward  the 
numerical  analysis  of  finite  open  structures.  The  advent  of 
the  method  of  moments  [14]  stimulated  the  frequency-domain 
study  of  these  structures;  while  recent  interest  in  tran¬ 
sient  methods  was  primarily  stimulated  by  a  study  due  to 
Bennett,  et  al  [43].  In  Bennett's  study,  both  transient 


99 


100 


numerical  and  experimental  results  for  several  canonic  open 
structures  were  presented.  Some  of  the  techniques  developed 
in  this  chapter  are  based  on  the  latter  contribution. 

The  study  of  the  SEM  parameters  for  open  structures  has 
been  limited  to  rectangular  geometries.  A  treatise  on  ob¬ 
taining  these  parameters  for  rectangular  apertures  and 
structures  using  frequency-domain  techniques  has  been  pre¬ 
sented  by  Pearson  [44]. 

In  this  chapter,  the  transient  numerical  solution  of  the 
thin,  perfect  conducting  rectangular  plate  problem  is  stu¬ 
died.  In  Section  6.2,  the  basic  mathematical  formulation  of 


the  problem  using  finite  difference  techniques  is  presented. 


The  basic  difference  formulation  is  then  act lied  to  a  'stan¬ 


dard'  gridding  scheme  in  Section  6.3,  and  a  'shifted'  grid- 
ding  scheme  in  Section  6.4.  The  shifted,  or  offset,  scheme 


was  initially  introduced  for  the  solution  of  rectangular 


problems  in  the  frequency-domain  by  Glisson  and  Wilton  [45]. 


Stability  analysis  and  current  distributions  for  each  of 


these  schemes  are  presented  in  the  appropriate  sections. 
Section  6.5,  the  TD-SEM  pole  distribution  for  the  square 
plate  is  introduced. 


6.2  MATHEMATICAL  FORMALISM 


The  electric  field  integral  expression  (2-12a)  is  the  ap¬ 
propriate  expression  for  describing  thin,  perfect  conducting 
structures.  For  the  geometry  shown  in  Figure  6.1,  this  ex¬ 
pression  may  be  written,  in  vector  potential  notation,  as 
the  following  coupled  set: 
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.  r  ,»x  »y  ,  =inc  ,  _inc  _inc  _inc  . 

where  A=(A  ,A7,0),  and  E  =(EX  ,E  , Ez  ). 
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The  following  representation  of  c.  will  be  used  for  all 


results  presented  in  this  chapter: 


EinC  (r;t)  *  0  E  exp  (~g2( (t-t)  +  k-r/c)2) 
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(6-2) 


where 


6  x  cosdcosii  +  y  cos8sin<|)  -  z  sin0,  (6-2a) 

k  «  x  3in0cos<p  +  y  3in0siiuj>  +  z  cos6. 


(6-2b) 


The  value  of  Ea  was  chosen  to  be  the  free  space  impedance. 


120ir,  t_„  was  chosen  to  be  2.45  light  meters  (LM),  and  g 
max 

-inc 

was  chosen  so  that  the  magnitude  of  E  at  t=0  and  t=4.9  LM 
would  be  0.0001  volts/meter. 

From  expansion  (3-17),  the  x  and  y  components  of  the  vec¬ 
tor  potential  may  be  explicitly  represented  as 
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and  A^’yn =AX’7  ( mA ,  nA ;  p A t ) . 

Unfortunately,  the  integral  which  appears  may  not  be 
evaluated  exactly  due  to  the  presence  of  the  time  pulse. 
Numerical  integration  or  linear  interpolation  represent  the 
possible  methods  of  evaluation.  In  Figure  6.2,  the  physical 
interpretation  of  how  the  time  pulse  activates  various  annu¬ 
lar  regions  which  contribute  to  the  integral  is  shown  (cau¬ 
sality  allows  us  to  only  consider  zero  or  positive  values  of 
the  time  difference  p-p' ) .  From  this  figure,  the  following 
linear  interpolation  formula  may  be  derived: 
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Int(D)  denotes  'integer  part  of',  and 
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3y  letting  a=(a*l/2),  b=(a-l/2),  c=( 3+1/2),  and  d=(3-l/2) 

we  may  evaluate  this  integral  as 
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A  simple  rectangular  rule  approximation  may  alternatively 


be  used  to  evaluate  integral  (6-5)  over  all  patches  other 
than  the  self  patch  without  introducing  considerable  error, 
i  .e. , 
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Use  of  linear  interpolation  and  approximation  (6-6)  is 
the  recommended  means  for  the  evaluation  of  equation  (6-3a). 
Numerical  integration  or  use  of  equation  (6-5a)  have  been 
found  to  introduce  only  slight  amplitude  shifts  in  the  final 
results,  and  therefore  their  use  is  unjustified  unless  pre¬ 
cise  results  are  sought.  The  interpretation  of  precise  is 
ambiguous,  however,  due  to  the  vast  number  of  models  which 
may  be  applied  to  a  particular  problem,  and  hence  the  vast 
number  of  slightly  different  results  which  may  be  obtained. 
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6.3  STANDARD  GRIPPING  SCHEME 

Figure  6.3  depicts  a  standard  gridding  scheme  for  a  rec¬ 
tangular  plate.  The  patches  as  shown  are  square  due  to  the 
choice  of  a  uniform  sampling  distance  in  both  the  x  and  y 
directions;  this  is  typically,  but  not  necessarily,  done. 

Two  complications  are  associated  with  this  simple,  commonly 
used,  model.  First,  since  the  two  current  density  compo¬ 
nents  lie  directly  on  top  of  one  another,  a  smooth  transi¬ 
tion  between  components  is  not  possible;  and  second,  boun¬ 
dary  conditions  on  the  current  density  must  be  explicitly 
enforced. 

For  the  standard  grid,  the  following  explicit  finite  dif¬ 
ference  scheme  for  the  x  component  of  the  vector  potential 
is  appropriate: 


AX  -  r2(A*  +  AX  ,  -  2AX  ) 

m.n.p+l  m-rl,n,p  m-l,n,p  m ,  n ,  p 


+  —  ( A^  +  A^  —  A^  —  \ 

4  m+l,n+l,p  m-l,n-l,p  m+l,n+l,p  m-l,n+l,p; 


+  A 


x 

®  *  Q  »  P  ™ 


1 


2A 


®  t  ^  9  P 


Einc 

X 


fn*0 , 1 
jn-0,1 
lP-1.2 


i 


(6-7) 


where  r2=(cAt/A)2.  Ay  may  be  similarly  defined.  An  expli¬ 
cit  difference  scheme  for  the  current  density  is  obtained  by 
substituting  expansion  (6-3)  into  the  above  difference  equa¬ 
tion  and  manipulating  the  indices  of  the  summations. 
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An  approximate  stability  criterion  for  the  vector  poten¬ 
tial  scheme  may  be  obtained  by  the  Fourier  stability  method 
(Sections  3.3.1,  5.3).  Let 


"ax" 

-  — 

a 

X 

X 

a 

L  yj 

•  « 

(»)  exp  •H(amA  +  8nA) 


(6-8) 


where  a  ,  a  are  arbitrary  coefficients,  u=exp{sAtj  (s  is 
x  y 


arbitrary  and  generally  complex),  and  a,  6  are  arbitrary 
real  constants.  By  substituting  the  expression  for  Ax  into 
equation  (6-7)  and  the  expression  for  Ay  into  a  similar 
equation,  we  obtain 


'  2 

OJ  *  u 


2-4r2sin2(aA/2)  r2sin(ctA)sin(SA)  1  f  a 


2  . 


2  0 

r'“sin(o(A)sin{sA)  2-4r  sir.4-  (  2A/2) 


J  l  yJ 


L 


-1  0 
0  -1 


1  |  a  . 

J  L  yj 


(6-9) 


This  matrix  equation  may  be  reduced  to  the  following  quartic 
equation  in  u: 


Ui4  +  (E+C)w3  +  (EC-BD+2)uj2  +  (E+C)w  +  1-0 


(6-10) 


where 


E  -  4r2sin2(aA/2)  -  2 


(6-10a) 


r-sin(aA)sin(2A) 


(6-10b) 


C  -  4r  sin  (BA/2) 


(6-10c) 


D  *  r^sin(0A)sin(aA)  . 


(6-10d) 


For  any  choice  of  real  a,  8,  |w|£l  for  rSl.  Therefore,  the 

choice  cAt=A  may  lead  to  a  stable  solution.  Certainty  can¬ 
not  be  obtained  since  boundary  conditions  cannot  be  included 
in  the  analysis.  The  matrix  stability  method  must  be  used 
to  verify  this  criterion  for  the  current  density  coefficient 
difference  formulation  with  arbitrary  boundary  conditions. 

The  enforcement  of  boundary  conditions  is  an  integral 
part  of  using  the  standard  gridding  scheme  effectively. 
Figure  6.4  shows  the  current  density  component  corresponding 
to  the  direction  of  polarization  of  the  incident  wave  on  the 
center  patch  of  a  one  meter  square  plate  discretized .with 
four  unknowns  in  each  direction  when  no  boundary  conditions 
are  enforced  and  choosing  cAt=0.7A.  The  result  is  highly 
oscillatory  yet  stable.  Figure  6.5  depicts  the  effect  of 
enforcing  components  of  the  current  perpendicular  to  the 
edges  to  be  zero.  In  Figure  6.6  an  attempt  has  been  made  to 


enforce  the  form  of  the  singular  behavior  of  the  current 


components  parallel  to  an  edge.  This  was  accomplished  by 
assuming  a  reciprocal  square  root  of  distance  singularity  as 
the  edge  is  approached  [42].  Hence,  the  value  /3  was  chosen 
as  the  extrapolation  constant  for  the  parallel  current  com¬ 
ponents  3A/2  from  the  edge,  i.e.,  the  parallel  component  at 
position  A/2  from  the  edge  is  V3  times  the  parallel  current 
component  evaluated  at  3A/^  from  the  edge. 

In  Figure  6.6  a  comparison  has  been  made  with  a  current 
density  distribution  obtained  from  the  theoretical  model 
used  by  3ennett  [43].  3ennett's  model  and  the  model  used 
here  differ  only  in  the  extrapolation  constant  used  for  the 
current  density  component  parallel  to  an  edge.  In  the  model 
used  by  3ennett,  an  extrapolation  constant  of  3  was  used  in¬ 
stead  of  /3 .  The  factor  of  3  was  found  to  occasionally 
yield  unstable  results,  whereas  /3  was  found  to  always  yield 
stable  results.  Hence,  the  latter  was  preferred.  The  two 
curves  agree  quite  closely  within  the  twelve  light  meter 
frame  which  is  shown.  In  should  be  noted  that  3ennett's  mo¬ 
del  has  been  shown  [43]  to  yield  results  quite  similar  to 
experimental  measurements. 

The  use  of  either  3  or  /3  as  an  extrapolation  technique 
i3  somewhat  unsatisfying  since  it  does  not  permit  a  time 
fluctuation  of  che  particular  form  which  is  being  forced. 

The  half  plane  problem  which  can  be  solved  analytically 


yields  a  reciprocal  square  root  behavior  for  the  current 
density  which  includes  both  a  spatial  and  a  time  dependence 
under  the  root  [42].  We  would  anticipate  a  similar  depen¬ 
dence  for  the  plate,  and  hence  any  extrapolation  should  re¬ 
flect  this.  Further  study  is  required. 

Although  theoretical  results  similar  to  experimental  mea 
surements  may  be  obtained  from  a  standard  gridding  scheme, 
it  is  difficult  to  know  a  priori  when  satisfactory  results 
have  been  obtained  due  to  the  modifications  which  are  re¬ 
quired.  In  an  effort  to  avoid  these  modifications  (i.e., 
create  a  more  natural  model),  we  consider  an  offset  or 
shifted  gridding  scheme. 

6.4  SHIFTED  GRIDDING  SCHEME 

Figure  6.7  depicts  a  shifted  gridding  scheme.  Three  de¬ 
sirable  properties  about  this  formulation  are  as  follows: 
current  is  allowed  to  make  a  smooth  transition  between  com¬ 
ponents,  zero  boundary  conditions  on  the  current  are  impli¬ 
citly  enforced,  and  it  is  not  necessary  to  step  off  the 
structure  for  any  finite  difference  evaluations. 

The  explicit  difference  representation  of  the  x  componen 
of  the  vector  potential  is  given  by 
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(6-11) 


A  similar  equation  may  be  developed  for  Ay .  An  explicit 
scheme  for  the  current  density  coefficients  may  be  obtained 
by  substituting  equation  (6-3)  into  the  above  equation  and 
manipulating  the  indices  of  the  summations. 

A  preliminary  stability  criterion  for  the  vector  poten¬ 
tial  scheme  may  be  found  by  the  Fourier  stability  method. 

By  letting 


s 

a  <jP  exp  {j(cmA  t  3(a-is)il)} 
x  I 

1 

a^P  exp  {j(a(m-is)A  +  Sni)  } 

(6-12) 


and  substituting  into  the  difference  scheme  (6-11),  we  ob¬ 
tain  the  following  quartic  equation  in  u  (similar  to  the  un¬ 
shifted  development) : 
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a4  +  (E+C)a)3  +  (EC-BEH-2)uj2  +  (E+C)u  +  1-0  (6-13) 

where 

E  -  4r2sin2(aA/2)  -  2 

B  -  r2(l+e-ja4ejSA-e-jaA-ej3A) 

(6-13a) 

C  -  4r2sin2(BA/2)  -  2 

D  -  r2(l+e-J6AejaA-e'j0A-ejaA). 

For  any  choice  of  real  a,  3,  l«|£l  for  r2S(l/2)  or 

cAt£(/2/2)A.  This  result  may  be  established  analyically  by 
assuming  four  solutions  of  the  form  expftjBj},  exp[±j02|  and 
noting  that  the  product  of  all  the  roots  must  be  unity,  or 
it  may  be  established  numerically. 

The  basic  result  may  be  extended  to  accommodate  three- 
dimensional  problems  with  different  spatial  sampling  distances 
in  each  of  the  three  spatial  directions,  i.e.,  we  have  Ax, 

Ay,  Az  instead  of  simply  A.  The  stability  requirement  for 
shifted  schemes,  in  general,  is  then 

(cAt)2  —  (Ax)2  +  (Ay)2  +  (Az)2  •  (6-14) 

The  validity  of  these  expressions  for  the  analogous  cur¬ 
rent  density  difference  formulation  must  be  confirmed  by  the 
matrix  stability  method. 


Figure  6.8  shows  the  unstable  result  obtained  by  choosing 
cAt=A.  By  choosing  cAt£(i/2/2)A,  however,  we  obtain  the  sta¬ 
ble  curves  shown  in  Figure  6.9.  A  comparision  of  these 
curves  with  the  curve  generated  by  the  standard  gridding 
scheme  has  been  made.  Note  that  the  amplitude  of  the  curves 
obtained  from  the  two  schemes  differ  slightly.  This  is  pri¬ 
marily  due  to  the  different  techniques  used  to  enforce  the 
zero  boundary  condition  in  each  scheme. 

6.5  TD-SEM  POLE  DISTRIBUTION 

In  this  section,  we  present  pole  distributions  obtained 
by  TD-SEM  (Section  3.3.1)  using  the  shifted  gridding  scheme 
on  a  one  meter  square  plate.  Figure  6.10  shows  the  distri¬ 
bution  for  a  total  of  2  unknowns  (A=l/2)  for  each  current 
component.  Figure  6.11  shows  the  distribution  for  a  total  of 
6  unknowns  (A=l/3),  and  Figure  6.12  shows  the  distribution 
for  a  total  of  12  unknowns  (A=l/4).  The  choice  cAt=0.7A  was 
made  throughout. 

As  wa3  discussed  in  Section  5. 3. 2.1,  an  ambiguity  exists 
in  the  validity  of  ail  the  'poles'  TD-SEM  yields  when  cAt  is 
chosen  less  than  A  due  to  an  increase  in  the  order  of  the 


transition  matrix.  A  filtering  scheme  to  remove  poles  which 
were  conjectured  to  be  a  consequence  of  the  numerical  pro- 
ceedure  was  discussed  in  that  section.  For  Figures  6.10-12 
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gure  6.11:  Poles  for  a  1  meter  square  plate  discretized 
with  a  shifted  gridding  scheme.  The  total 
number  of  patches  for  each  component  of  the 
current  density  was  6,  and  the  choice  cAt=0.7A 
was  made. 
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Figure  6.12:  Poles  for  a  1  meter  square  piate  discretized 
with  a  shifted  gridding  scheme.  The  total 
number  of  patches  for  each  component  of  the 
current  density  was  12,  and  the  choice  cAt=0.7A 
was  made. 


a  similar,  but  slightly  different,  filtering  scheme  was 
used.  The  scaled  poles  corresponding  to  the  eigenvalues 
with  positive  real  components  were  restricted  to  a  maximum 
magnitude  on  the  imaginary  axis  of  1/(2A)  radians;  the  sca¬ 
led  poles  corresponding  to  the  eigenvalues  with  negative 
real  components  were  restricted  to  the  range  1/(2A)  to  1/A 
radians.  A  line  has  been  drawn  on  these  figures  to  separate 
the  two  regions  since  the  true  physical  meaning  of  the  poles 
corresponding  to  eigenvalues  with  negative  real  components 
is  not  clear  for  rectangular  geometries  (although  it  was 
conjectured  that  these  poles  have  no  meaning  for  the  wire 
problem) .  A  sensible  method  to  test  the  validity  of  the 
poles  in  both  regions  is  to  reduce  the  value  of  cAt  below 
the  initial  choice  of  0.7A  and  note  shifts  in  the  pole  posi¬ 
tions.  For  true  poles,  we  suspect  very  little  shift.  It 
was  experimentally  observed  that  the  lower  set  shifted  only 
slightly,  but  the  upper  set  experienced  a  considerable 
shift.  Hence  from  this  argument,  we  conjecture  that  only 
the  lower  set  represents  true  system  poles. 

The  lower  pole  cluster  agrees  reasonably  well  with  the 
frequency-domain  results.  The  lowest  order  pole  from  Figure 
6.12  (12  unknowns)  is  explicitly  -0.284+j0.715;  this  pole 
was  found  by  frequency-domain  methods  to  be  -0 . 272- jO . 675 . 
The  time-domain  result  should  approach  this  (up  to  the  simi¬ 
larity  of  the  models  used  in  each  domain)  as  the  number  of 


unknowns  is  increased.  Figure  6.13  compares  selected 
time-domain  poles  with  available  frequency-domain  poles. 

It  is  interesting  to  note  that  double  poles  were  observed 
for  the  square  plate,  and  a  fill  in  of  poles  occured  as  the 
number  of  unknowns  increased.  Both  of  these  can  be  justi¬ 
fied  by  considering  the  poles  of  an  infinite  rectangular 
waveguide . 

The  simultaneous  iteration  method  presented  in  Chapter  4 
may  be  used  to  obtain  the  natural  frequencies  and  modes  for 
the  rectangular  plate.  Care  is  required  in  its  implementa¬ 
tion,  however,  due  to  the  absence  of  the  well  defined  layer¬ 
ing  structure  which  appeared  with  the  wire  problem. 


Chapter  VII 
CONCLUSIONS 

In  this  thesis,  the  fundamental  integral  equations  of 
electromagnetic  theory  and  the  theoretical  foundations  of 
both  the  frequency-domain  SEM  and  the  time-domain  SEM  method 
of  Cordaro  and  Davis  were  developed.  It  was  observed  from 
the  development  of  TD-SEM  that  special  sparse  eigensolution 
routines  were  required  to  determine  the  eigenvalues  of  the 
transition  matrix  due  to  the  excessive  high  speed  storage 
requirements  associated  with  certain  problems.  A  modified 
simultaneous  iteration  algorithm  was  developed  to  satisfy 
this  eigensolution  requirement.  The  algorithm  may  be  used 
to  obtain  partial  pole  solutions  for  a  variety  of  geome¬ 
tries,  and  was  explicitly  shown  to  be  effective  on  the 
thin-wire  problem  for  an  arbitrary  number  of  unknowns. 

Root  searching  methods  which  take  advantage  of  the  com¬ 
panion  form  of  the  transition  matrix,  such  as  Muller's  meth¬ 
od  and  the  polynomial  matrix  reduction  method  of  Woolivich, 
were  found  to  be  effective  methods  for  obtaining  the  natural 
frequencies  only  for  linear  geometries  discretized  with  re¬ 
latively  few  unknowns.  The  contour  integration  technique  of 
Singaraju,  Giri,  and  Baum,  which  also  takes  advantage  of 
this  companion  form,  was  not  explicitly  tested  in  this  stu- 
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dy;  however,  the  method  may  prove  to  be  effective  for  prob¬ 
lems  discretized  with  a  large  number  of  unknowns  provided 
the  required  integration  operation  can  benefically  be  used 
to  numerically  smooth  the  determinate  evaluation.  Further 
study  is  required  to  establish  the  utility  of  this  method. 

Time-domain  techniques  which  provide  convenient  matrix 
methods  for  obtaining  the  SEM  coupling  coefficients  have 
been  developed.  The  time-domain  form  of  these  coefficients 
is  much  simpler  than  the  equivalent  frequency- domain  form. 

The  effect:  of  altering  the  sub-sectional  coupling  between 
unknowns  in  the  numerical  formulation  of  the  thin-wire  prob¬ 
lem  was  also  investigated.  No  relation  between  the  specific 
layering  structure  of  the  poles  and  sub-matrices  of  the 
transition  matrix  was  observed.  This  was  unfortunate  since 
it  was  hoped  that  if  only  a  particular  subset  of  the  pole 
distribution  was  desired,  then  sub-matrices  which  did  not 
influence  this  subset  could  be  removed  from  the  transition 
matrix  and  thereby  yield  a  lower  order  problem. 

A  shifted  gridding  scheme  was  applied  to  the  rectangular 
plate  problem  to  obtain  transient  solutions.  This  scheme 
was  found  to  represent  a  more  natural  discretization  for  the 
problem  than  the  unshifted  or  standard  gridding  scheme  which 
is  typically  used.  The  shifted  scheme  was  then  used  in  con¬ 
junction  with  TD-SEM  to  obtain  pole  distributions  for  the 
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square,  perfect  conducting  plate.  The  results  obtained  were 
found  to  be  quite  similar  to  available  frequency-domain  re¬ 
sults. 

Fourier  and  matrix  stability  methods  were  applied  to  the 
finite  difference  representations  of  electromagnetic  equa¬ 
tions.  These  methods  were  found  to  provide  accurate  insight 
into  the  required  relation  between  the  time  and  spatial  sam¬ 
pling  distances  that  will  yield  a  numerically  stable  solu¬ 
tion  for  an  arbitrary  difference  formulation. 

The  physical  significance  of  the  additional  poles  gener¬ 
ated  by  choosing  the  time  sampling  distance  smaller  than  the 
spatial  sampling  distance  or  the  spatial  sampling  distance 
to  be  different  in  different  directions  remains  an  open 
question.  These  additional  poles  are  conjectured  to  be 
false  poles  for  which  an  elimination  procedure  has  been  pre¬ 


sented. 
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